c Version 2.51 from April 30, 1998 	
c********************************************************************* 
      SUBROUTINE ISMCO(PARM,OBJ,PHCON,PHREF,PHBLB,NE,IIGEOM,IIDISTF,
     + IIREFL) 
c********************************************************************* 
c this program calculates the photon flux 
c for the Compton scattered component (PHCON, normalized to 1 at 1 keV) 
c Reflected component (PHREF), and Black body (PHBLB)  
c for the disk-corona model neglecting energy balance 
c for the given parameters of the problem: 
c A: Maxwellian electrons 
c  TE (keV)  - electron temperature of the corona
c B: Power-law electrons
c  PLIND     - index of the power-law 
c  GMIN      - minimum Lorentz factor  
c  GMAX      - maximum Lorentz factor  
c C: Maxwellian + power-law tail
c  TE (keV)  - electron temperature of the corona
c  PLIND     - index of the power-law
c  GMIN      - Lorentz factor   of the break 
c  GMAX      - maximum Lorentz factor  
c 
c  TPH (keV) - temperature of the initial radiation (black-body) 
c  TAU       - vertical optical depth of  corona 
c  ETA       - cosine of the inclination angle (if <0 rescattering) 
c  GCOVER    - covering factor or Rin inner radius of the disc 
c  REFLEC    - amount of reflection 
c  FEAB      - iron abundance in units of solar
c  REDSHIFT  - redshift 
c  ab_       - abundance of elements heavier than He relative to
c              the solar abundances
c  xi        - disk ionization parameter = L/nR^2
c  temp      - disk temperature (for ionization balance, Tmax=10^6 K) 
c  beta      - reflection emissivity law (r^beta)
c  radin     - inner radius of the disk
c  radout    - outer radius of the disk
c  t0taug    - ratio of the height to the radius of the cylinder 
c********************************************************************* 
      IMPLICIT REAL*8(A-H,O-Z)
      logical flag,flagfre
      integer NE,IIGEOM,IIDISTF,IIREFL
      PARAMETER(IDEVC=12) 
      PARAMETER(NDATA=10) 
C NG - number of points in Gauss quadrature 
C (used in   POWMATR & MAXGALA, MAXCUT)  
      PARAMETER(NG=5)
      PARAMETER(NS=5)
C NGO - number of points in Gauss quadrature 
C (used in OUTFLUX) 
      PARAMETER(NGO=3)
C NC=MAXANG -  number of points in Gauss quadrature (used in 
C          angular integrals) 
C NT=MAXTAU - number of points in trapec formula 
C            (integration over  optical depth)   
      PARAMETER(MAXTAU=15)
      PARAMETER(nem=5000) 
      PARAMETER(np=20,nrefl=7)     
      PARAMETER(MAXFRE=71,MAXANG=3,MAXSC=50)
      PARAMETER(II=MAXFRE,NSP=2*II,NC=MAXANG,ND=1,KK=II*NC,LL=ND*KK,
     1NLYR=MAXTAU,NL1=NLYR+1)
      PARAMETER(NEX=II+II*(II-1)/2,NCX=NC+NC*(NC-1)/2)
      PARAMETER(NX1=II+1)
      PARAMETER(IPC=ND*II)
      COMMON/JPWFRE/XX(II),XLOG(II),XKEV(II),EXXY(NEX),S0(II)
      COMMON/QQX12/XX12(NX1) 
      COMMON/QWITCH/IGROUND,ICHAND,KEYRAY,IREDF,IDISTF,IGEOM,IPOL,ISP  
      COMMON/SWITC2/EPSMIN,XEMAX,NISO,ISOTR
      COMMON/JPKN/TPH,TE,THOM,HX,HPHK02,NL,NT,NSCAT,NPOL  
      COMMON/QREFLEC/UREF(LL),DINTEN(II),TOPESC(KK),SIDESC(KK) 
      COMMON/QQRES/DINTPL(MAXSC,LL),DINTMI(MAXSC,LL), 
     *SUIPL(LL),SUIMI(LL),SUIMIOUT(LL) 
      COMMON/QQRST/FRDPP(LL,LL),FRDPM(LL,LL),FRDX(IPC,IPC)
      COMMON/QQCM2/A(II),UANG(NC),AANG(NC),AC(LL),AINT(LL),CINT(LL)
      COMMON/QQWG/UG(NG),AG(NG)
      COMMON/QQWGO/UGO(NGO),AGO(NGO) 
      COMMON/QQWS/UGS(NS),AGS(NS)
      COMMON/GEOM/ETAIN,ETAOUT,ETACOV 
      COMMON/QQLUM/GSOFT,GCOVER,FILF1,FILFAC 
      COMMON/QWPH/YPH,YDISK 
      COMMON/QPOWE/XMIN,XMAX,POW  
      COMMON /JPPOW/PLIND,GMIN,GMAX,DNORPOW,DNORMAX,RATIO 
      COMMON /QQWYY/ YE,Y1,YR1,YR2,DY,YKY,DK02,DK2Y,G0 

      DIMENSION DSPXX(II),F2(II),SUIDAVE(II),OBSANG(1),OBSREF(NEM)  
      DIMENSION PARM(np),FINT(NC),pa0(np)
      real*4 par_refl(nrefl)
      DIMENSION OBSCON(NEM,NC),OBSBLB(NEM)
      DIMENSION PHCON(NEM),PHREF(NEM),PHBLB(NEM),OBJ(NEM)
      DIMENSION BBINT(NEM),XFRE(II),OBJOLD(NEM) 
      DIMENSION DIMULT(NEM),SUIDOWN(LL),IOB(NEM)
      save    pa0,flag,objold,neold,flagfre
      save    compabs,obscon,obsblb,obsref,bbint
      save    igeold,iob,phnorm,nkev,ionold 
      data    pa0/np*9999./
      data    flag/.true./,flagfre/.true./
      data    neold/0/
      DATA CKEV/0.0019569472D0/,PI/3.141592653589792384D0/,EPS/1D-14/
      DATA C13/0.3333333333333333D0/,C38/0.375D0/,C316PI/5.9683104D-2/ 
      DATA CLGMIN/-30d0/,CMIN/1d-30/
      DATA IN/5/,NOUT/6/,EPSREF/1D-6/  
C ************************************
C ************************************

      IF(IIDISTF.eq.0) THEN 
      TE=DABS(PARM(1)/511.)
      ELSEIF(IIDISTF.eq.2) THEN 
      PLIND=DABS(PARM(2))
      GMIN=PARM(3)
      GMAX=PARM(4)
      ELSEIF(IIDISTF.eq.4) THEN
      TE=DABS(PARM(1)/511.)
      PLIND=PARM(2) 
      GMIN=PARM(3)
      GMAX=PARM(4) 
      ENDIF 
      TPH=DABS(PARM(5)) 
      TAU=DABS(PARM(6)) 
      ETA=DABS(PARM(7))
      IF(IIGEOM.le.2) THEN
      GCOVER=PARM(8)
      ELSE 
      RIN=PARM(8)
      ENDIF 
      REFLEC=PARM(9)
      FEAB=DABS(PARM(10)) 
      REDSHIFT=PARM(14)
      if(IIREFL.le.1) then
      IONIZ=IIREFL 
      IREFL=0
      par_refl(2)=0. 
      else 
      IONIZ=IIREFL-2  
c relativistic disk reflector 
      IREFL=1
c  ab_       - abundance of elements heavier than He relative to
c              the solar abundances
      par_refl(2)=alog10(sngl(PARM(11)))
      endif 
c reflection parameters 
c Fe abudance 
      par_refl(1)=alog10(sngl(FEAB)) 
c  xi        - disk ionization parameter = L/nR^2
      par_refl(3)=sngl(PARM(12))
c  temp      - disk temperature (for ionization balance, Tmax=10^6 K) 
      par_refl(4)=sngl(PARM(13))
c      temp=sngl(TPH)*1.1604448e7 
c  beta      - reflection emissivity law (r^beta)
      par_refl(5)=sngl(PARM(15))
c  radin     - inner radius of the disk
      par_refl(6)=sngl(PARM(16))
c  radout    - outer radius of the disk
      par_refl(7)=sngl(PARM(17))
C IGEOM=0 SLAB geometry
C IGEOM=1 CYLINDER geometry (with T0TAUG=H/R)
C IGEOM=2 HEMISPHERE geometry  (TOTAUG=1, H=R) 
c for slab  GIN=1, for ARs: cylinder H=R: g<0.6; H=2R: g<0.5; 
c for hemisphere g<0.7
       if(IIGEOM.le.2) THEN 
       IGEOM=IIGEOM 
       IGROUND=0
       GRUNC=1D0-GCOVER 
       FILFAC=1D0
       FILF1=0D0
       else 
       IGEOM=IIGEOM-3
       IGROUND=2
       GCOVER=1. 
       IF(IGEOM.eq.0) GIN=1.
       IF(IGEOM.eq.1) GIN=0.5
       IF(IGEOM.eq.2) GIN=0.7
       IF(IGEOM.eq.0) GOUT=0. 
       IF(IGEOM.eq.1) GOUT=0.1
       IF(IGEOM.eq.2) GOUT=0.05 
       PTAU=1D0-DEXP(-TAU) 
       FILF1=RIN*RIN 
       IF(IGEOM.eq.0) FILF1=1D0-RIN 
       FILFAC=1D0-FILF1 
       endif 

C********************************
c     if the energy vector changed from the previous call
      if(ne.ne.neold)then
         flagfre=.true.
         neold=ne
      else 
         do 10 i=1,ne
            if(obj(i).ne.objold(i)) then
               flagfre=.true.
            endif 
 10      continue 
      endif

      if(flagfre) then 
         do 11 i=1,ne
            objold(i)=obj(i)
 11      continue
      endif 

C************************************
C INPUT PARAMETERS 
C************************************
c EPSMIN accuracy of the spectrum computed 
c XEMAX  maximum energy up to what accuracy is checked 
      EPSMIN=3.D-3
      XEMAX=2.D0 
c ICOLD=0 normal coronal radiation 
      ICOLD=0 
c      IF(PARM(7).lt.0D0) ICOLD=1 

c T0TAUG=H/R 
c for slab T0TAUG=0; for cylinder =parm(18),  for hemisphere =1
      IF(IGEOM.eq.0) T0TAUG=0.
      IF(IGEOM.eq.1) T0TAUG=PARM(18) 
      IF(IGEOM.eq.2) T0TAUG=1.
C IREDF =0 exact redistribution function    
C IREDF =1 isotropic approximation in electron rest frame 
      IF(PARM(1).gt.0..or.PARM(2).gt.0.) THEN 
      IREDF=0
      ELSE 
      IREDF=1
      ENDIF  
C NISO=0 - exact scattering process
C NISO=1 - anisotropy for first 0,..ISOTR-1 scatterings,
C          for all other, isotropic TAU dependent source function
      IF(PARM(6).gt.0.) THEN 
      NISO=1
      ISOTR=3
      ELSE 
      NISO=1
      ISOTR=1 
      ENDIF  
C XLMIN - decimal logarithms of low  photon energy 
C               cutoff (in units of electron rest mass)  
c XLMAX - maximum of the spectral range 
c 0.5 eV 
      XLMIN=-6.0D0
c 50 MeV 
      XLMAX=2.0D0
C*****************************************
C END OF MAIN INPUT PARAMETERS 
C*****************************************
      MAXITER=1 
      IF(REFLEC.lt.EPSREF) MAXITER=0 
      ISP=0 
C IDISTF=0 Maxwellian electron distribution 
C IDISTF=1 monoenergetic electron distribution 
C IDISTF=2 power-law  electron distribution 
C IDISTF=4 Maxwellian + power-law tail
      IDISTF=IIDISTF 
C***************************************
C
C other parameters: 
C ICHAND=0 isotropic distribution of initial radiation 
C ICHAND=1 Chandrasekhar-Sobolev distribution of intensity (1+2*mu)
C ICHAND=2 Isotropic flux I=1/mu 
      ICHAND=0 
C NL    -  number of points in Gauss-Laguerre quadrature  
      NL=10
      NT=MAXTAU
      NX=II
C
C *****************************************       
C
C constants 
      DLN=DLOG(1D1)  
      C38PI=C38/(PI*PI) 
      PI2=2D0*PI
C------------------------------
      YPH=511D0/TPH
      ETA2=ETA*ETA
      WETA=1D0-ETA
      WETA2=WETA*WETA

C************************************

      if(flag) then 
      CALL AGUEGALA(NL,EPS) 
      CALL QDRGSDO(UANG,AANG,NC)  
      DO 20 L=1,NC
         UANG(L)=0.5D0*(1D0+UANG(L))                         
         AANG(L)=AANG(L)*5D-1                                
   20 CONTINUE 
      CALL QDRGSDO(UG,AG,NG)  
      DO 23 L=1,NG
         UG(L)=0.5D0*(1D0+UG(L))                         
         AG(L)=AG(L)*5D-1                                
   23 CONTINUE 
      CALL QDRGSDO(UGO,AGO,NGO)  
      DO 25 L=1,NGO
         UGO(L)=0.5D0*(1D0+UGO(L))                         
         AGO(L)=AGO(L)*5D-1                                
   25 CONTINUE 
      CALL QDRGSDO(UGS,AGS,NS)
      DO 26 L=1,NS
         UGS(L)=0.5D0*(1D0+UGS(L))
         AGS(L)=AGS(L)*5D-1
   26 CONTINUE
c frequency step 
         HX=(XLMAX-XLMIN)/DFLOAT(NX-1)
         CONS=DLN*HX 
C calculation of frequencies
         CALL FREQPOIN(XLMIN,XLMAX)
C weights of frequency integration 
      DO 27 JX=1,II 
      AJ=1D0
      IF(JX .eq. 1 .or. JX .eq. II) AJ=5D-1 
      A(JX)=AJ*CONS
   27 CONTINUE
c weights of frequency-angular integration 
      DO 29 N=1,ND
      DO 29 I=1,II
      DO 29 J=1,NC
      K=J+(I-1)*NC 
      L=(N-1)*KK+K
c AINT and CINT are FREQUENCY and ANGLE  supervectors.
      AINT(L)=A(I)
      CINT(L)=AANG(J)
      AC(L)=A(I)*AANG(J) 
   29 CONTINUE                                           
      flag=.false.
      endif 

      if(pa0(14).ne.parm(14)) then 
c XFRE - vector of energies in the observer's frame 
      DO 195 IX=1,II  
      XFRE(IX)=XX(IX)/(1D0+REDSHIFT) 
  195 continue 
c  CKEV=1 keV in units of mc2  
      CALL DLOCATE(XFRE,II,CKEV,NKEV)
      endif 

C blackbody intrinsic radiation (isotropic incident radiation) 
      if(pa0(5).ne.parm(5)) then
      if(PARM(5).ge.0d0) CALL BLACKB(XX,DINTEN,II) 
      if(PARM(5).lt.0d0) CALL BBMULTI(XX,DINTEN,II) 
      endif 

C******************************************************
C******************************************************
C******************************************************
c initial values for reflected spectrum 
      DO 37 L=1,LL 
      UREF(L)=0D0 
   37 CONTINUE 
      IF(IDISTF.eq.0.or.IDISTF.eq.4) YE=1D0/TE 
      if(pa0(1).ne.parm(1).or.pa0(2).ne.parm(2).or.
     & pa0(3).ne.parm(3).or.pa0(4).ne.parm(4))then
      CALL COMWY(EPS) 
C Compton cross-section 
      CALL CROSCOMP 
c COMPTON scattering redistribution matrix 
      CALL CSRF 
      endif 

      if(flagfre) then 
      DO 41 I=1,NE 
      CALL DLOCATE(XX,II,OBJ(I),IOB(I))
   41 CONTINUE 
      endif 

      NUMBITER=0  
 
      GSOFT=GCOVER
  
       if(pa0(1).eq.parm(1).and.pa0(2).eq.parm(2).and.
     &  pa0(3).eq.parm(3).and.pa0(4).eq.parm(4)
     &  .and.pa0(5).eq.parm(5).and.pa0(6).eq.parm(6).and.
     &  pa0(9).eq.parm(9).and.pa0(10).eq.parm(10).and. 
     &  pa0(8).eq.parm(8).and.igeom.eq.igeold.and.
     &  ioniz.eq.ionold.and.pa0(13).eq.parm(13).and.
     &  pa0(11).eq.parm(11).and.pa0(12).eq.parm(12).and.
     &  pa0(14).eq.parm(14).and.pa0(15).eq.parm(15).and. 
     &  pa0(16).eq.parm(16).and.pa0(17).eq.parm(17).and. 
     &  pa0(18).eq.parm(18)) then
      if(flagfre) then 
cjp energies changed 
      goto 2000
      else
cjp only angle changed 
      goto 3000 
      endif 
      endif 
C******************************************************
       ETAIN=GIN/2D0
       ETAOUT=0.5D0-ETAIN
       ETACOV=ETAIN*FILFAC+ETAOUT
C beginning of iteration procedure 
   50 CONTINUE         
C******************************************************
C calculation of the comptonized spectrum 
      CALL COMPSPEC(TAU,T0TAUG)
C 
      DO 54 L=1,LL
      if(IGEOM.ne.0) then 
      SUIDOWN(L)=FILFAC*SUIMI(L)+FILF1*SUIMIOUT(L)
      else 
      SUIDOWN(L)=SUIMI(L) 
      endif 
   54 CONTINUE
c flux down averaged over angles
      DO 212 I=1,II
      sum=0d0
      DO 210 IANG=1,NC
      K=IANG+(I-1)*NC
      sum=sum+aang(iang)*SUIDOWN(K)*uang(iang)
  210 CONTINUE
      SUIDAVE(I)=sum 
  212 CONTINUE
c 
      if(REFLEC.ge.EPSREF) then
C  CALCULATION OF THE REFLECTED COMPONENT 
      if(ioniz.eq.0.and.PARM(10).lt.0D0) then 
c Poutanen, Nagendra & Svensson 
      call PNSrefl(suidown,xx,uref,uang,
     + par_refl,nrefl,ii,kk,nc,ioniz,0)
      else 
c Magdziarz & Zdziarski 1995 
c using  angle averaged incoming flux 
      call MZrefl(xx,suidave,xx,uref,uang,
     + par_refl,nrefl,ii,ii,nc,kk,ioniz,0)
      endif 
      endif 

      SUCM=0D0
      SUCP=0D0 
      SUCMOUT=0D0 
      SUREF=0D0 
      DO 135 IX=1,NX 
      X=XX(IX)
cjp300697 
      AJ=A(IX) 
      AX=AJ*X 
c      AX=A(IX)
c      AJ=AX/X 
      SUMCM=0D0
      SUMCP=0D0
      SUMCMOUT=0D0 
      SUMREF=0D0 
      DO 137 IANG=1,NC 
      K=IANG+(IX-1)*NC 
      UZ=UANG(IANG)
      AZ=AANG(IANG)
c comptonized bb+ref 
      SUMCM=SUMCM+AZ*SUIMI(K)*UZ 
      SUMCP=SUMCP+AZ*SUIPL(K)*UZ 
      SUMCMOUT=SUMCMOUT+AZ*SUIMIOUT(K)*UZ
c initial reflected 
      SUMREF=SUMREF+AZ*UREF(K)*UZ  
  137 CONTINUE  
      SUREF=SUREF+SUMREF*AX
      SUCM=SUCM+SUMCM*AX
      SUCP=SUCP+SUMCP*AX 
      SUCMOUT=SUCMOUT+SUMCMOUT*AX 
  135 CONTINUE 
      DLC=SUCP+SUCM 
      AETA=SUCM/DLC
      ETAOUT=SUCMOUT/DLC
      ETAIN=AETA-ETAOUT
      ETACOV=ETAIN*FILFAC+ETAOUT 
      SUCM=PI2*SUCM
      SUREF=PI2*SUREF
      IF(IGEOM.eq.1.or.IGEOM.eq.2) THEN
      SUCM=PI*SUCM
      SUREF=PI*SUREF
      ENDIF
c calculation of the absorption compactness 
      COMPABS=SUCM-SUREF
c multiplication with REFLEC
cjp no multiplication with REFLEC for the reflection which is comptonized
       DO 130 L=1,LL 
       if(IGEOM.eq.0) UREF(L)=UREF(L)*FILFAC
c       if(IGEOM.ne.0) UREF(L)=UREF(L)*REFLEC
  130  CONTINUE 

      NUMBITER=NUMBITER+1  
      IF(NUMBITER .GT. MAXITER) GO TO 140
      GO TO 50  
C****************************************************
C end of iteration procedure 
C ****************************************************
  140 continue
c**********************************
 2000 CONTINUE  
C END of ITERATION PROCEDURE 
C ********************************************

C ****************************************************
C
      if(ICOLD.eq.0) then
c calculation of COMPTON SCATTERED FLUX at obs. grid 
      DO 150 IANG=1,NC 
      DO 152 I=1,II 
      K=IANG+(I-1)*NC 
      DSPXX(I)=DLOG(SUIPL(K)+CMIN)  
  152 CONTINUE 
      CALL QSPLINE(XLOG,DSPXX,II,1d30,1d30,F2) 
      DO 150 I=1,NE 
      IF(IOB(I).lt.1.or.IOB(I).ge.II) THEN 
      OBSCON(I,IANG)=CMIN 
      ELSE 
      CALL QSPLINT(XLOG,DSPXX,F2,II,DLOG10(OBJ(I)),YYY)
      OBSCON(I,IANG)=DEXP(YYY) 
      ENDIF 
  150 CONTINUE 

c calculation of BLACK BODY RADIATION at  obs. grid 
      if(PARM(5).ge.0d0) CALL BLACKB(OBJ,BBINT,NE) 
      if(PARM(5).lt.0d0) CALL BBMULTI(OBJ,BBINT,NE)
c multi-colour disk at  obs. grid
      CALL BBMULTI(OBJ,DIMULT,NE) 

c reflection at obs. grid
C using Magdziarz & Zdziarski 1995
      if(REFLEC.lt.EPSREF) then 
      DO 205 I=1,NE
      OBSREF(I)=CMIN
  205 CONTINUE 
      else 
c neutral reflector 
      obsang(1)=ETA

      if(ioniz.eq.0.and.PARM(10).lt.0D0) then 
c Poutanen, Nagendra & Svensson 
      call PNSrefl(suidown,obj,obsref,obsang,
     + par_refl,nrefl,ne,ne,1,ioniz,0)
      else 
c Magdziarz & Zdziarski 1995 
c using  angle averaged incoming flux 
      call MZrefl(xx,suidave,obj,obsref,obsang,
     + par_refl,nrefl,ii,ne,1,ne,ioniz,irefl)
      endif 

c multiplication with REFLEC 
      do I=1,NE 
      if(IGEOM.ne.0) OBSREF(I)=OBSREF(I)*REFLEC
      enddo 
C end of reflection  at obs. grid 
      endif 

c TAKING INTO ACCOUNT COVERING FACTOR 

      DO 170 IX=1,NE
      DO 172 IANG=1,NC
      K=IANG+(IOB(IX)-1)*NC 
      IF(IOB(IX).lt.1.or.IOB(IX).ge.II) THEN 
       IF(IOB(IX).eq.0) K=IANG  
       ESCPROB=TOPESC(K)+SIDESC(K) 
      ELSE 
       K1=IANG+IOB(IX)*NC 
       ESCMIN=TOPESC(K)+SIDESC(K)
       ESCMAX=TOPESC(K1)+SIDESC(K1)
       CALL DLINEAR(XX(IOB(IX)),ESCMIN,XX(IOB(IX)+1),ESCMAX,
     * OBJ(IX),ESCPROB)  
      ENDIF 
      FINT(IANG)=ESCPROB
  172 CONTINUE 
      CALL POLINTQ(UANG,FINT,NC,ETA,FLY,DELY)
      ESCPROB=FLY
      IF(IGROUND.eq.0) THEN  
      OBSBLB(IX)=(GRUNC*COMPABS+GSOFT*ESCPROB)*BBINT(IX)*
     * REFLEC  
      OBSREF(IX)=(GRUNC+GSOFT*ESCPROB)* 
     * OBSREF(IX) 
      ELSE 
        if(IGEOM.eq.1.or.IGEOM.eq.2) then 
      OBSBLB(IX)=DIMULT(IX)*(1D0-PTAU*GOUT)+
     +               ESCPROB*BBINT(IX)*FILFAC  
      OBSREF(IX)=(ETAOUT+ETAIN*FILFAC*ESCPROB)/ETACOV*
     *               OBSREF(IX) 
        else 
      OBSBLB(IX)=ESCPROB*BBINT(IX)
      OBSREF(IX)=(FILFAC*ESCPROB+
     + DMAX1(0D0,REFLEC-FILFAC))*OBSREF(IX) 
        endif 
      ENDIF 
  170 CONTINUE 
      endif 

 3000 continue 
C normalization of the Compton scattered 
c photon spectrum  (=1 at observer's 1 keV) 
c interpolation of SUIPL to the observed angle
      IF(ICOLD.eq.0) THEN 
      DO 204 NK=NKEV,NKEV+1
      DO 202 IANG=1,NC
      K=IANG+(NK-1)*NC 
      FINT(IANG)=SUIPL(K)*UANG(IANG)
  202 continue
      CALL POLINTQ(UANG,FINT,NC,ETA,FLY,DELY)
      IF(NK.eq.NKEV)DEL1=FLY/XX(NKEV) 
      IF(NK.eq.NKEV+1)DEL2=FLY/XX(NKEV+1) 
  204 continue
      CALL DLINEAR(XFRE(NKEV),DEL1,XFRE(NKEV+1),DEL2, 
     *            CKEV,PHNORM) 

c interpolation between angles to get photon flux at ETA 
      DO 190 IX=1,NE  
      DO 182 IANG=1,NC
      FINT(IANG)=OBSCON(IX,IANG)*UANG(IANG)  
  182 continue 
      CALL POLINTQ(UANG,FINT,NC,ETA,FLY,DELY)
      PHCON(IX)=FLY/OBJ(IX) 
      PHBLB(IX)=OBSBLB(IX)/OBJ(IX)*ETA 
      PHREF(IX)=OBSREF(IX)/OBJ(IX)*ETA  
  190 continue 
      ENDIF 


c renormalization to 1 at 1 keV 
      DO 200 IX=1,NE  
      PHCON(IX)=PHCON(IX)/PHNORM 
      PHBLB(IX)=PHBLB(IX)/PHNORM 
      PHREF(IX)=PHREF(IX)/PHNORM 
  200 continue 
C
      flagfre=.false. 
      ionold=ioniz
      igeold=igeom 

         do n=1,np
            pa0(n)=parm(n)
         enddo
      RETURN 
      END                  
c***********************************************************************
c        REDISTRIBUTION MATRICES FOR COMPTON SCATTERING 
c***********************************************************************
c***********************************************************************
      SUBROUTINE CSRF 
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER (IDEVC=12)
      PARAMETER(NG=5) 
      PARAMETER(NS=5) 
      PARAMETER(MAXTAU=15) 
      PARAMETER(MAXFRE=71,MAXANG=3,MAXSC=50)
      PARAMETER(II=MAXFRE,NSP=2*II,NC=MAXANG,ND=1,KK=II*NC,LL=ND*KK,
     1NLYR=MAXTAU,NL1=NLYR+1)
      PARAMETER(NEX=II+II*(II-1)/2,NCX=NC+NC*(NC-1)/2)
      PARAMETER(IPC=ND*II)
      COMMON/QQCM2/A(II),UANG(NC),AANG(NC),AC(LL),AINT(LL),CINT(LL)
      COMMON/QQRST/FRDPP(LL,LL),FRDPM(LL,LL),FRDX(IPC,IPC)
      COMMON/JPKN/TPH,TE,THOM,HX,HPHK02,NL,NT,NSCAT,NPOL 
      COMMON /QQWYY/ YE,Y1,YR1,YR2,DY,YKY,DK02,DK2Y,G0 
      COMMON/QWITCH/IGROUND,ICHAND,KEYRAY,IREDF,IDISTF,IGEOM,IPOL,ISP 
      COMMON/JPWFRE/XX(II),XLOG(II),XKEV(II),EXXY(NEX),S0(II)
      COMMON/QQWG/UG(NG),AG(NG) 
      COMMON/QQWS/UGS(NS),AGS(NS) 
      DIMENSION COMPHI11(II),REDXX(NS),DMUA(NS),F2(NS)   
      DATA PI/3.141592653589792384D0/,EPS/1D-14/ 
      DATA C13/0.33333333333333333D0/,NOUT/6/ 
      NX=II 
c-------------------------------------------------------------------------
C CALCULATIONS OF ANGLE DEPENDENT REDISTRIBUTION FUNCTIONS EXPLICITELY FOR
c (CSM+RT) CODE. 
c-------------------------------------------------------------------------
      DLN=DLOG(1D1) 
      H3=PI*2D0 
      DO 11 IAU=1,NS 
      DMUA(IAU)=UGS(IAU)*2D0-1D0 
   11 CONTINUE 
c------------------------------------------------------------------
c begin  STAGE I (portions taken from subroutine PHIMATR in csm2.f)
C mainly integrals over azimuth angle:
C mainly integrals over azimuth angle:
      DO 1020 L1=1,LL
      DO 1020 L=1,LL
      FRDPP(L,L1)=0.D0
      FRDPM(L,L1)=0.D0
1020  CONTINUE
      JXMIN=1
        DO 400 IX=1,NX
        X=XX(IX)
          IF(IDISTF.eq.0) JXMIN=IX
          DO 300 JX=JXMIN,NX
          X1=XX(JX)
          CALL COMX(X,X1)
c here computation of R(x,x1,mu) 
             DO 350 IAU=1,NS 
             QMU=1D0-DMUA(IAU) 
             IF(QMU.gt.EPS) THEN 
             CALL REDFUNC(X,X1,QMU,SMU) 
             ELSE 
             SMU=0D0
             ENDIF 
             REDXX(IAU)=SMU  
  350 CONTINUE 
             CALL QSPLINE(DMUA,REDXX,NS,1d30,1d30,F2)

            DO 200 KETA=1,NC  
            K=KETA+(IX-1)*NC   
            ETA=UANG(KETA) 
            ETA2=ETA*ETA
            SI2=1D0-ETA2
            SI=DSQRT(SI2) 
              DO 100 LETA=1,KETA
c              DO 100 LETA=1,NC 
              K1=LETA+(JX-1)*NC   
              ETA1=UANG(LETA)
              ETA12=ETA1*ETA1
              SI12=1D0-ETA12
              SI1=DSQRT(SI12)  
              SIKL=SI*SI1
              COKL=ETA*ETA1
              SUM11P=0D0 
              SUM11M=0D0
c BEGIN DELTA (azimuth) INTEGRATION LOOP
c                 PHI=0D0
c                 NSIGN=-1
                 DO 50 ID=1,NS
                 PHI=UGS(ID)*PI 
c                 NA=NSIGN+3 
c                 IF(ID.eq.1.or.ID.eq.NS) NA=1 
                 CFI=DCOS(PHI)
                 PHKL=CFI*SIKL
c BEGIN CALN. OF  S+- (+mu,-mu') = S-+ (-mu,+mu') BLOCK MATRIX:
                 DMUPM=-COKL+PHKL
                 CALL QSPLINT(DMUA,REDXX,F2,NS,DMUPM,SM)
c                 CALL REDFUNC(X,X1,QMUPM,SM) 
c END CALN. OF  S+- (+mu,-mu') = S-+ (-mu,+mu') BLOCK MATRIX.
c BEGIN CALN. OF  S++ (+mu,+mu') = S-- (-mu,-mu') BLOCK MATRIX:
                 DMUPP=COKL+PHKL
                 CALL QSPLINT(DMUA,REDXX,F2,NS,DMUPP,SP)
c                 CALL REDFUNC(X,X1,QMUPP,SP) 
c END CALN. OF  S++ (+mu,+mu') = S-- (-mu,-mu') BLOCK MATRIX.
                 SUM11P=SUM11P+SP*AGS(ID)
                 SUM11M=SUM11M+SM*AGS(ID)
c                 NSIGN=-NSIGN
c                 PHI=PHI+HPHI
   50 CONTINUE              
c END  DELTA (azimuth) INTEGRATION LOOP. 
c 
c--------------------------------------------------------------------------
      FRDPP(K1,K)=DMAX1(H3*SUM11P,0D0) 
      FRDPM(K1,K)=DMAX1(H3*SUM11M,0D0)
c CONSTRUCTION OF A PART OF THE REDISTRIBUTION MATRICES FRD's IS OVER.
  100 CONTINUE 
  200 CONTINUE 
  300 CONTINUE 
  400 CONTINUE 
c end of STAGE I
c---------------
c------------------------------------------------------------------
c begin  STAGE II 
c CONSTRUCTION OF FULL STRUCTURE OF FRD's IS DONE BELOW BY EXPLOITING THE
c SYMMETRY PROPERTIES OF THE CSM.
c ATTENTION! DO NOT CHANGE THE ORDER of LOOPS !!! 
c power-law electrons 
      IF(IDISTF.NE.0) THEN
      DO 180 I=1,NX
      X=XX(I)
      DO 170 J=1,NC
      K=J+(I-1)*NC
      DO 160 I1=1,NX
      X1=XX(I1)
      XXR=X/X1
      KE1=J+(I1-1)*NC
      DO 150 J1=1,NC
      KE=J1+(I-1)*NC
      K1=J1+(I1-1)*NC
        IF(J1.gt.J) THEN 
        FRDPP(K1,K)=FRDPP(KE1,KE)*XXR
        FRDPM(K1,K)=FRDPM(KE1,KE)*XXR
        ELSE
        FRDPP(K1,K)=FRDPP(K1,K)*XXR
        FRDPM(K1,K)=FRDPM(K1,K)*XXR
        ENDIF
  150 CONTINUE
c END OF mu' LOOP
  160 CONTINUE
c END OF X' LOOP
  170 CONTINUE
c END OF mu LOOP
  180 CONTINUE
c END OF X LOOP
c--------------------------------------------------------------------------
      ELSE 
c Maxwellian electrons 
      DO 280 I1=1,NX
      X1=XX(I1)
      DO 270 I=1,NX
      X=XX(I)
      XXR=X/X1
      JSX=I1+I*(I-1)/2
      IF(I.GT.I1)  EXXJS=EXXY(JSX)*XXR
      DO 260 J=1,NC
      K=J+(I-1)*NC
      KE1=J+(I1-1)*NC
      DO 250 J1=1,NC
      KE=J1+(I-1)*NC
      K1=J1+(I1-1)*NC
        IF(I.LE.I1) THEN
         IF(J1.gt.J) THEN
         FRDPP(K1,K)=FRDPP(KE1,KE)*XXR
         FRDPM(K1,K)=FRDPM(KE1,KE)*XXR
         ELSE
         FRDPP(K1,K)=FRDPP(K1,K)*XXR
         FRDPM(K1,K)=FRDPM(K1,K)*XXR
         ENDIF
        ELSE 
         IF(J1.lt.J) THEN
         FRDPP(K1,K)=FRDPP(KE,KE1)*EXXJS
         FRDPM(K1,K)=FRDPM(KE,KE1)*EXXJS
         ELSE 
         FRDPP(K1,K)=FRDPP(K,K1)*EXXJS
         FRDPM(K1,K)=FRDPM(K,K1)*EXXJS
         ENDIF 
        ENDIF 
c--------------------------------------------------------------------------
c--------------------------------------------------------------------------
  250 CONTINUE
c END OF mu' LOOP
  260 CONTINUE
c END OF mu LOOP
  270 CONTINUE
c END OF X LOOP
  280 CONTINUE
c END OF X' LOOP
      ENDIF 
c NOW WE SPECIFICALLY USE THE SYMMETRY RELATIONS: FRDMP=FRDPM; FRDMM=FRDPP:
c end of STAGE II
c begin  STAGE III
c CHECKING THE NORMALIZATION OF THE ANGLE DEPENDENT GENERAL  CSM 
c REDISTRIBUTION MATRIX:
c FIRST DEVELOP THE ANGLE AVERAGED REDISTRIBUTION MATRICES FOR WORK 
      DO 835 I=1,NX
      DO 835 I1=1,NX
      sum=0d0 
      DO 834 J=1,NC
      K=J+(I-1)*NC
      DO 834 J1=1,NC
      K1=J1+(I1-1)*NC
c      IF(L.eq.L1) FRDPP(L,L1)=FRDPP(L,L1) + SCON(J)
      sum=sum+ CINT(K)*(FRDPP(K1,K)+FRDPM(K1,K))*CINT(K1)
834   CONTINUE
      FRDX(I1,I)=sum 
835   CONTINUE
c FRDX GIVEN ABOVE IS AN ANGLE AVERAGED REDISTRIBUTION FUNCTION. 
      DO 840 I=1,NX
      X=XX(I)
      SUM=0.D0
      DO 841 I1=1,NX
      AX1=XX(I1)*A(I1) 
      SUM=SUM + FRDX(I,I1)*AX1
841   CONTINUE
      COMPHI11(I)=SUM 
840   CONTINUE
c-----------------------------------------------------------------------------
c COMPHI(NX) GIVEN ABOVE IS <computed CS profile function>. IN PRINCIPLE,
c IT SHOULD BE SAME AS  S0(NX)  EVALUATED DIRECTLY IN THE CSM PART OF THE CODE
c (KN CROSS SECTION BASICALLY). BUT IN PRACTICE IT IS <not>, DUE TO THE FINITE
c QUADRATURE APPROXIMATIONS IN PERFORMING THE NUMERICAL INTEGRATIONS. HENCE 
c THERE IS A NEED TO RENORMALIZE THE REDISTRIBUTION FUNCTION.
c-----------------------------------------------------------------------------
      DO 849 I=1,NX 
      DO 849 I1=1,NX
      CONS=S0(I1)/COMPHI11(I1)
c NOW THE RENORMALIZATION OF THE ANGLE AVERAGED 
      FRDX(I1,I)=FRDX(I1,I) * CONS
849   CONTINUE
      DO 850 I=1,NX 
      DO 850 J=1,NC 
      K=J+(I-1)*NC
      DO 850 I1=1,NX
      CONS=S0(I1)/COMPHI11(I1)
      DO 850 J1=1,NC
      K1=J1+(I1-1)*NC
c NOW THE RENORMALIZATION OF THE FULL ANGLE DEPENDENT  <CSM>
c REDISTRIBUTION MATRICES  FRDPP,FRDPM 
      FRDPP(K1,K)=FRDPP(K1,K) * CONS
      FRDPM(K1,K)=FRDPM(K1,K) * CONS
850   CONTINUE
c END  OF THE OF CALCULATION OF NORMALIZED CSM REDISTRIBUTION MATRICES 
c [vector S_{ij} (++)] = [vector S_{ij} (--)]   AND
c [vector S_{ij} (+-)] = [vector S_{ij} (-+)],  i,j = 1,2.
c WE ARE NOW READY TO EMPLOY THESE REDISTRIBUTION MATRICES IN (CSM+RT).
c end of STAGE III
c----------------
      RETURN
      END
C**********************************************************************
C
C IGEOM =  0 SLAB     
C IGEOM =  1 CYLINDER 
C IGEOM =  2 HEMISHPERE 
C integration over the optical depth
C iteration prosedure:  Neuman series. 
C 
      SUBROUTINE COMPSPEC(TAU,T0TAUG)
      IMPLICIT REAL*8(A-H,O-Z)           
      PARAMETER(IDEVC=12)
      PARAMETER(MAXTAU=15)
      PARAMETER(MAXFRE=71,MAXANG=3,MAXSC=50)
      PARAMETER(II=MAXFRE,NSP=2*II,NC=MAXANG,ND=1,KK=II*NC,LL=ND*KK,
     1NLYR=MAXTAU,NL1=NLYR+1)
      PARAMETER(NEX=II+II*(II-1)/2,NCX=NC+NC*(NC-1)/2)
      PARAMETER(NX1=II+1)
      PARAMETER(IPC=ND*II)
      COMMON/JPWFRE/XX(II),XLOG(II),XKEV(II),EXXY(NEX),S0(II)
      COMMON/QWITCH/IGROUND,ICHAND,KEYRAY,IREDF,IDISTF,IGEOM,IPOL,ISP
      COMMON/SWITC2/EPSMIN,XEMAX,NISO,ISOTR
      COMMON/JPKN/TPH,TE,THOM,HX,HPHK02,NL,NT,NSCAT,NPOL 
      COMMON/QWPH/YPH,YDISK 
      COMMON/QQX12/XX12(NX1) 
      COMMON/GEOM/ETAIN,ETAOUT,ETACOV 
      COMMON/QREFLEC/UREF(LL),DINTEN(II),TOPESC(KK),SIDESC(KK)   
      COMMON/QQLUM/GSOFT,GCOVER,FILF1,FILFAC
      COMMON/QQRST/FRDPP(LL,LL),FRDPM(LL,LL),FRDX(IPC,IPC)
      COMMON/QQRES/DINTPL(MAXSC,LL),DINTMI(MAXSC,LL), 
     *SUIPL(LL),SUIMI(LL),SUIMIOUT(LL) 
      COMMON/QQCM2/A(II),UANG(NC),AANG(NC),AC(LL),AINT(LL),CINT(LL)
      COMMON /QQWYY/ Y,YR,YR1,YR2,DY,YKY,DK02,DK2Y,G0 
      COMMON/WAGUE/UL(40),AL(40),AU(40),AE(40),N0        
      DIMENSION SOUPL(MAXTAU,LL),SOUMI(MAXTAU,LL), 
     1 BGPL(MAXTAU,MAXTAU,KK),BGMI(MAXTAU,MAXTAU,KK), 
     2 BGPL1(MAXTAU,MAXTAU,KK),BGMI1(MAXTAU,MAXTAU,KK),   
     3 BGPL2(MAXTAU,MAXTAU,KK),BGMI2(MAXTAU,MAXTAU,KK), 
     4 BGMIN(MAXTAU,MAXTAU,KK),BGMEX(MAXTAU,MAXTAU,KK), 
     5 BGPLA(MAXTAU,MAXTAU,II),BGMIA(MAXTAU,MAXTAU,II),
     6 BGMINA(MAXTAU,MAXTAU,II)
      DIMENSION TAUVEC(MAXTAU),SOUPL0(MAXTAU,LL),SOUMI0(MAXTAU,LL), 
     1 DITAPL0(LL,MAXTAU),DITAMI0(LL,MAXTAU)  
      DIMENSION DOUT(MAXTAU,KK),SOISTA(MAXTAU,II),DIMETA(II,MAXTAU)
      DIMENSION CINP(MAXTAU,MAXTAU,NC),CINM(MAXTAU,MAXTAU,NC), 
     1          CEMP(MAXTAU,MAXTAU,NC),CEMM(MAXTAU,MAXTAU,NC)
      save TAUOLD,YOLD,FIL1OLD,IGEOLD,YPHOLD,DOUT  
      save BGPL,BGMI,BGPL1,BGMI1,BGPL2,BGMI2,BGMIN,BGMEX
      save BGPLA,BGMIA,BGMINA 
     
      DATA CKEV/0.0019569472D0/,PI/3.141592653589792384D0/,EPS/1D-14/ 
      DATA C13/0.33333333333333333D0/,C12/5D-1/
      DATA NOUT/6/ 
      data TAUOLD/9D9/,YOLD/-9D9/,FIL1OLD/9D9/,IGEOLD/10/,YPHOLD/-9D9/  
      DO 377 K=1,KK 
      SUIPL(K)=0D0 
      SUIMI(K)=0D0 
 377  CONTINUE 

      PI2=2D0/PI 
      T0TA2=T0TAUG/2D0
      DNT=DFLOAT(NT-1)
      DNTR=1D0/DNT 
      DELTAU=TAU*DNTR 
      DLN=DLOG(1D1)

      DTT=0D0
      DO 4 ITAU=1,NT 
      TAUVEC(ITAU)=DTT
      DTT=DTT+DNTR
    4 CONTINUE 

      if(IGEOM.ne.IGEOLD.or.YOLD.ne.Y.or.TAUOLD.ne.TAU) then 
c AUXILIARY MASSIVES 
      DO 3 IX=1,II   
      DO 3 IANG=1,NC 
      K=IANG+(IX-1)*NC
      ETA=UANG(IANG) 
      SIET=DSQRT(1D0-ETA*ETA) 
      TSI=PI2*T0TAUG*SIET 
      DO 3 ITAU=1,NT 
      BGPL(ITAU,ITAU,K)=1D0 
      BGMI(ITAU,ITAU,K)=1D0 
      BGPL1(ITAU,ITAU,K)=1D0 
      BGMI1(ITAU,ITAU,K)=1D0 
      IF(IGEOM.eq.1) THEN 
      BGPL2(ITAU,ITAU,K)=TSI 
      BGMI2(ITAU,ITAU,K)=TSI 
      ENDIF 
      IF(IGEOM.eq.2) THEN 
      DTR=TAUVEC(ITAU)
      ETADTR=ETA*DTR 
      BGPL2(ITAU,ITAU,K)=ETADTR 
      IF(DTR.gt.SIET)  BGPL2(ITAU,ITAU,K)=2D0*ETADTR 
      BGMI2(ITAU,ITAU,K)=ETADTR
      ENDIF 
    3 CONTINUE 

      DO 9 IANG=1,NC 
      ETA=UANG(IANG) 
      SIET=DSQRT(1D0-ETA*ETA) 
      SIETA=SIET/ETA
      TAUETA=TAU/ETA
      TSI=PI2*T0TAUG*SIET 
      DO 7 ITAU=2,NT 
      ITAU1=ITAU-1 
      DTI=TAUVEC(ITAU)
      ETADTI=ETA*DTI
      DO 5 LTAU=1,ITAU1 
      DTL=TAUVEC(LTAU)
      CINPL=0D0 
      CINMI=0D0 
      CEMPL=0D0 
      CEMMI=0D0 
      TT0=DABS(DTL-DTI)*SIETA 
      IF(IGEOM.eq.0 .or. IGEOM.eq.1) THEN 
      TTAG=TT0*T0TA2 
      IF(TTAG .lt. 1D0) THEN 
      CEMPL=DSQRT(1D0-TTAG**2) 
      CINPL=(DACOS(TTAG)-TTAG*CEMPL)*PI2 
      CEMPL=CEMPL*TSI
      CINMI=CINPL
      CEMMI=CEMPL
      ENDIF 
      ENDIF 
      IF(IGEOM.eq.2) THEN 
      TT2=TT0*TT0
      IF(ITAU.eq.NT) THEN 
      RI=0D0  
      RI2=0D0 
      ELSE 
      RI2=1D0-DTI*DTI  
      RI=DSQRT(RI2)
      ENDIF
      RL2=1D0-DTL*DTL  
      RL=DSQRT(RL2)
      IF(TT0.le.(RL-RI)) THEN 
      CINPL=1D0 
      CINMI=RI2/RL2
      CEMPL=2D0*ETADTI 
      CEMMI=0D0 
      ENDIF 
      IF(TT0.gt.(RL-RI).and.TT0.lt.(RI+RL)) THEN 
      TT0RI2=2D0*TT0*RI
      TT0RL2=2D0*TT0*RL
      PHIST=DACOS((RI2-RL2+TT2)/TT0RI2) 
      PHIDU=DACOS((RL2-RI2+TT2)/TT0RL2) 
      DQRR=DSQRT(4D0*RI2*RL2-(TT2-RI2-RL2)**2) 
      CPL=(RI2*PHIST+RL2*PHIDU-0.5D0*DQRR)/PI
      CINPL=CPL/RI2 
      CINMI=CPL/RL2    
      SINSTI=DQRR/TT0RI2 
c      IF(ETA.ge.ETATANI) THEN 
      PHI0=PHIST
      SINPHI0=SINSTI 
c      ELSE 
c      PHI0=DMIN1(PHIST,PI/2D0) 
c      SINPHI0=DMIN1(1D0,SINSTI) 
c      ENDIF 
      CEMPL=PI2*(ETADTI*PHI0+SIET*RI*SINPHI0)
      SINSTL=DQRR/TT0RL2 
      CEMMI=PI2*(-ETA*DTL*PHIDU+SIET*RL*SINSTL)
      ENDIF 
      ENDIF 
      CINP(LTAU,ITAU,IANG)=CINPL
      CINM(LTAU,ITAU,IANG)=CINMI
      CEMP(LTAU,ITAU,IANG)=CEMPL
      CEMM(LTAU,ITAU,IANG)=CEMMI
    5 CONTINUE 
    7 CONTINUE 
    9 CONTINUE 


      DO 19 IX=1,II 
      S0IX=S0(IX) 
      DO 17 IANG=1,NC 
      K=IANG+(IX-1)*NC 
      ETA=UANG(IANG) 
      TAUETA=TAU/ETA
      DO 15 ITAU=2,NT 
      ITAU1=ITAU-1 
      DTI=TAUVEC(ITAU)
      ETADTI=ETA*DTI
      DO 13 LTAU=1,ITAU1 
      DTL=TAUVEC(LTAU)
      TDDE=TAUETA*(DTI-DTL) 
      EGPL=0D0
      SUPL=S0IX*TDDE 
      IF(SUPL.lt.7d1) EGPL=DEXP(-SUPL)  
      BGPL(LTAU,ITAU,K)=EGPL
      BGMI(LTAU,ITAU,K)=EGPL
c for radiation field inside 
      BGPL1(LTAU,ITAU,K)=EGPL*CINP(LTAU,ITAU,IANG) 
      BGMI1(LTAU,ITAU,K)=EGPL*CINM(LTAU,ITAU,IANG) 
c for the emergent radiation field 
      BGPL2(LTAU,ITAU,K)=EGPL*CEMP(LTAU,ITAU,IANG) 
      BGMI2(LTAU,ITAU,K)=EGPL*CEMM(LTAU,ITAU,IANG) 
   13 CONTINUE 
   15 CONTINUE 
   17 CONTINUE 
   19 CONTINUE 
c angle averaged
      DO 21 ITAU=1,NT
      DO 21 LTAU=1,ITAU
      DO 21 IX=1,II
      SUM1=0D0
      SUM2=0D0
      DO 23 IANG=1,NC
      AZ=AANG(IANG)/UANG(IANG)
      K=IANG+(IX-1)*NC
      SUM1=SUM1+AZ*BGPL1(LTAU,ITAU,K)
      SUM2=SUM2+AZ*BGMI1(LTAU,ITAU,K)
   23 CONTINUE
      BGPLA(LTAU,ITAU,IX)=SUM1
      BGMIA(LTAU,ITAU,IX)=SUM2
   21 CONTINUE
      
c FOR THE DISK WITH HOLES
      IF(IGROUND.eq.2) THEN 
      DO 109 IANG=1,NC
      ETA=UANG(IANG)
      SIET=DSQRT(1D0-ETA*ETA)
      SIETA=SIET/ETA
      TAUETA=TAU/ETA
      DO 107 ITAU=1,NT
      DTI=TAUVEC(ITAU)
      ETADTI=ETA*DTI 
      DO 105 LTAU=1,NT
      DTL=TAUVEC(LTAU)
      TDDE=TAUETA*(DTI+DTL) 
      CINPL=0D0
      CEMPL=0D0
      TT0=(DTL+DTI)*SIETA

c SLAB AND CYLINDER
      IF(IGEOM.eq.0 .or. IGEOM.eq.1) THEN
      TTAG=TT0*T0TA2
      IF(TTAG .lt. 1D0) THEN
      CEMPL=DSQRT(1D0-TTAG**2)
      CINPL=(DACOS(TTAG)-TTAG*CEMPL)*PI2
      CEMPL=CEMPL*PI2*T0TAUG*SIET
      ENDIF
      ENDIF

c HEMISHPERE
      IF(IGEOM.eq.2) THEN
      TT2=TT0*TT0
      IF(ITAU.eq.NT) THEN
      RI2=0D0
      ELSE
      RI2=1D0-DTI*DTI
      ENDIF
      IF(LTAU.eq.NT) THEN
      RL2=0D0
      ELSE
      RL2=1D0-DTL*DTL
      ENDIF
      RI=DSQRT(RI2)
      RL=DSQRT(RL2)

      IF(TT0.le.DABS(RL-RI)) THEN
       IF(RL.ge.RI) THEN
      CINPL=1D0
      CEMPL=2D0*ETADTI
       ELSE
      CINPL=RL2/RI2
      CEMPL=0D0
       ENDIF
      ENDIF
      IF(TT0.gt.DABS(RL-RI).and.TT0.lt.(RI+RL)) THEN
      TT0RI2=2D0*TT0*RI
      TT0RL2=2D0*TT0*RL
      PHIST=DACOS((RI2-RL2+TT2)/TT0RI2) 
      PHIDU=DACOS((RL2-RI2+TT2)/TT0RL2) 
      DQRR=DSQRT(4D0*RI2*RL2-(TT2-RI2-RL2)**2)
      CPL=(RI2*PHIST+RL2*PHIDU-0.5D0*DQRR)/PI
      SINSTI=DQRR/TT0RI2
      CINPL=CPL/RI2
      CEMPL=PI2*(ETADTI*PHIST+SIET*RI*SINSTI)
      ENDIF
      ENDIF
      CINP(LTAU,ITAU,IANG)=CINPL
      CEMP(LTAU,ITAU,IANG)=CEMPL
  105 CONTINUE
  107 CONTINUE
  109 CONTINUE

      DO 119 IX=1,II
      S0IX=S0(IX) 
      DO 117 IANG=1,NC
      K=IANG+(IX-1)*NC
      ETA=UANG(IANG)
      TAUETA=TAU/ETA
      DO 115 ITAU=1,NT
      DTI=TAUVEC(ITAU)
      DO 113 LTAU=1,NT
      DTL=TAUVEC(LTAU)
      TDDE=TAUETA*(DTI+DTL) 
      CINPL=0D0
      SUPL=S0IX*TDDE 
      EGPL=0D0
      IF(SUPL.lt.7d1) EGPL=DEXP(-SUPL)
c for radiation field inside
      BGMIN(LTAU,ITAU,K)=EGPL*CINP(LTAU,ITAU,IANG)
c for the emergent radiation field
      BGMEX(LTAU,ITAU,K)=EGPL*CEMP(LTAU,ITAU,IANG)
  113 CONTINUE
  115 CONTINUE
  117 CONTINUE
  119 CONTINUE
c angle averaged 
      DO 121 ITAU=1,NT
      DO 121 LTAU=1,NT
      DO 121 IX=1,II 
      SUM1=0D0
      DO 123 IANG=1,NC 
      ETA=UANG(IANG)
      K=IANG+(IX-1)*NC
      SUM1=SUM1+AANG(IANG)/ETA*BGMIN(LTAU,ITAU,K) 
  123 CONTINUE 
      BGMINA(LTAU,ITAU,IX)=SUM1 
  121 CONTINUE 
      ENDIF

      DO 132 K=1,KK 
      IF(IGEOM.eq.2) THEN 
      TOPESC(K)=0D0
      ELSE 
      TOPESC(K)=BGPL1(1,NT,K) 
      ENDIF 
      IF(IGEOM.eq.0) SIDESC(K)=0D0 
  132 CONTINUE 

      IF(IGEOM.eq.1.or.IGEOM.eq.2) THEN    
c unscattered radiation emergent from sides 
        DO 110 IANG=1,NC
           ETA=UANG(IANG) 
           DNTETA=DNT*ETA 
        DO 110 IX=1,II
           K=IANG+(IX-1)*NC 

           SUPL0=0D0
c unscattered bb+ref
            DO 100 ITAU=1,NT 
            ATAU=1D0 
            IF(ITAU.eq.1.or.ITAU.eq.NT) ATAU=5D-1 
            SUPL0=SUPL0+ATAU*BGPL2(1,ITAU,K) 
  100 CONTINUE 
            SUPL0=SUPL0/DNTETA 
            SIDESC(K)=SUPL0 
  110 CONTINUE
      ENDIF

c end of TAU Y IGEOM condition 
      endif 

         ISC=1 
C outside soft input
      IF((IGROUND.eq.2).and.(IGEOM.eq.1.or.IGEOM.eq.2).and.
     & (IGEOM.ne.IGEOLD.or.TAUOLD.ne.TAU.or.YPHOLD.eq.YPH)) 
     & CALL OUTFLUX(TAUVEC,DOUT,TAU,T0TAUG)

C                                             
C UNSCATTERED RADIATION (HERE ISC=1)
C
C    LOOP ON THE FREQUENCY 
        DO 80 IX=1,II
        DINT=DINTEN(IX)
C    LOOP ON THE ANGLES
           DO 70 IANG=1,NC
           K=IANG+(IX-1)*NC 
            ETA=UANG(IANG) 
            EPTA=DELTAU/ETA
      IF(IGROUND.eq.0.or.IGEOM.eq.0) THEN
      DNAREF=DINT*QANGLE(ETA,ICHAND)+GCOVER*UREF(K)
      ELSE 
      DNAREF=FILFAC*(4D0*DINT*QANGLE(ETA,ICHAND)+ETAIN/ETACOV*UREF(K)) 
      ENDIF 
C    LOOP ON THE OPTICAL DEPTH 
              DO 60 ITAU=1,NT
              LTAU=NT-ITAU+1 
              DITAPL0(K,ITAU)=DNAREF*BGPL1(1,ITAU,K)
      IF(IGROUND.eq.2.and.(IGEOM.eq.1.or.IGEOM.eq.2)) 
     & DITAPL0(K,ITAU)=DITAPL0(K,ITAU)+DOUT(ITAU,K)
              DITAMI0(K,LTAU)=0D0 
      SOUPL(ITAU,K)=0D0 
      SOUMI(ITAU,K)=0D0
   60 CONTINUE 
      SUIMIOUT(K)=0D0 
   70 CONTINUE 
   80 CONTINUE 
  
      ISC=2 
      DIFMAX=1D0 
C MULTIPLY SCATTERING: ITERATION PROCEDURE 
      DO WHILE(ISC .lt. MAXSC .and. DIFMAX .gt. EPSMIN) 
      IF(ISC.gt.2) DIFMAX=0D0 

c calculation of the mean intensity dependent on TAU 
      IF((NISO.eq.1).and.(ISC.eq.ISOTR+1)) THEN 
        DO 415 ITAU=1,NT
        DO 415 IX=1,II
         SUMTI=0D0
        DO 419 IANG=1,NC
        K=IANG+(IX-1)*NC
        SUMTI=SUMTI+AANG(IANG)*(DITAPL0(K,ITAU)+DITAMI0(K,ITAU))
  419 CONTINUE
      DIMETA(IX,ITAU)=SUMTI*0.5D0
  415 CONTINUE
      ENDIF 

c ISOTROPIC SOURCE FUNCTIONS dependent on TAU
      IF(NISO.eq.1.and.ISC.ge.ISOTR+1) THEN
      DO 444 ITAU=1,NT
      DO 444 IX=1,II
      X=XX(IX)
      SUMI=0D0
      DO 445 JX=1,II
      SUMI=SUMI+FRDX(JX,IX)*DIMETA(JX,ITAU)*A(JX)
  445 CONTINUE
      SOISTA(ITAU,IX)=SUMI*X
  444 CONTINUE
      DO 446 IX=1,II
      X=XX(IX) 
      DO 446 IANG=1,NC
      L=IANG+(IX-1)*NC 
      DO 446 ITAU=1,NT
      XXPM=SOISTA(ITAU,IX)
      SOUPL0(ITAU,L)=XXPM
      SOUMI0(ITAU,L)=XXPM 
      SOUPL(ITAU,L)=XXPM+SOUPL(ITAU,L)
      SOUMI(ITAU,L)=XXPM+SOUMI(ITAU,L)
           IF(ITAU.eq.NT.and.XXPM.gt.1D-30.and.X.lt.XEMAX) THEN
           DIF=XXPM/SOUPL(NT,L)
           DIFMAX=DMAX1(DIF,DIFMAX)
           ENDIF
  446 CONTINUE
      ENDIF
c ACCURATE SOURCE FUNCTIONS
      IF(NISO.eq.0.or.ISC.lt.ISOTR+1) THEN 
        DO 190 IX=1,II 
        X=XX(IX) 
        DO 190 IANG=1,NC
        L=IANG+(IX-1)*NC 
        DO 185 ITAU=1,NT
           XIPL=0D0 
           XIMI=0D0 
        DO 180 JX=1,II 
        DO 180 JANG=1,NC
        L1=JANG+(JX-1)*NC 
        FRDPPEL=FRDPP(L1,L)
        FRDPMEL=FRDPM(L1,L)
        DIACPL=DITAPL0(L1,ITAU)*AC(L1)
        DIACMI=DITAMI0(L1,ITAU)*AC(L1)
        XIPL=XIPL+FRDPPEL*DIACPL+FRDPMEL*DIACMI
        XIMI=XIMI+FRDPMEL*DIACPL+FRDPPEL*DIACMI 
  180 CONTINUE 
            XXPL=XIPL*X 
            XXMI=XIMI*X  
            SOUPL0(ITAU,L)=XXPL
            SOUMI0(ITAU,L)=XXMI
            SOUPL(ITAU,L)=XXPL+SOUPL(ITAU,L) 
            SOUMI(ITAU,L)=XXMI+SOUMI(ITAU,L)  
  185 CONTINUE 
       IF(ISC.gt.2.and.XXPL.gt.1D-30.and.X.lt.XEMAX) THEN
       DIF=XXPL/SOUPL(NT,L)
       DIFMAX=DMAX1(DIF,DIFMAX)
       ENDIF
  190 CONTINUE 
      ENDIF

c INTERNAL RADIATION FIELD 
      IF(NISO.eq.0.or.ISC.lt.ISOTR+1) THEN 
        DO 280 IX=1,II 
        DO 280 IANG=1,NC
        K=IANG+(IX-1)*NC 
        L=K 
          ETA=UANG(IANG)
          EPTA=DELTAU/ETA 
              DO 260 ITAU=1,NT
              LTAU=NT-ITAU+1 
              SUPL1=0D0 
            IF(IGROUND.eq.2) THEN
                DO 251 IT=1,NT
                AT=1D0
                IF(IT.eq.1.or.IT.eq.NT) AT=5D-1
                SUPL1=SUPL1+AT*SOUMI0(IT,L)*BGMIN(IT,ITAU,K)
  251 CONTINUE
            SUPL1=SUPL1*FILF1 
            ENDIF
            IF(ITAU.ne.1) THEN 
                DO 250 IT=1,ITAU 
                AT=1D0 
                IF(IT.eq.1.or.IT.eq.ITAU) AT=5D-1 
                SUPL1=SUPL1+AT*SOUPL0(IT,L)*BGPL1(IT,ITAU,K)
  250 CONTINUE 
            ENDIF 
              SUPL1=SUPL1*EPTA 
              SUMI1=0D0 
            IF(LTAU.ne.NT) THEN 
                DO 240 IT=LTAU,NT 
                AT=1D0 
                IF(IT.eq.LTAU.or.IT.eq.NT) AT=5D-1 
                SUMI1=SUMI1+AT*SOUMI0(IT,L)*BGMI1(LTAU,IT,K)
  240 CONTINUE 
              SUMI1=SUMI1*EPTA 
            ENDIF
              DITAPL0(L,ITAU)=SUPL1
              DITAMI0(L,LTAU)=SUMI1 
  260 CONTINUE 
c emergent radiation field up-down 
      IF(IGEOM.ne.2) SUIPL(L)=SUIPL(L)+DITAPL0(L,NT)
                     SUIMI(L)=SUIMI(L)+DITAMI0(L,1)
  280 CONTINUE 
      ELSE 
c for angle averaged radiation field 
c emergent radiation field up-down
        DO 301 IX=1,II
        DO 303 IANG=1,NC
        K=IANG+(IX-1)*NC
        L=K
          ETA=UANG(IANG)
          EPTA=DELTAU/ETA
              SUPL1=0D0
              SUPL2=0D0
      IF(IGEOM.ne.2)THEN 
                DO 305 IT=1,NT
                AT=1D0
                IF(IT.eq.1.or.IT.eq.NT) AT=5D-1
      IF(IGROUND.eq.2) SUPL2=SUPL2+AT*SOUMI0(IT,L)*BGMIN(IT,NT,K)
                       SUPL1=SUPL1+AT*SOUPL0(IT,L)*BGPL1(IT,NT,K)
  305 CONTINUE
              SUPL1=(SUPL1+SUPL2*FILF1)*EPTA  
      ENDIF 
              SUMI1=0D0
                DO 307 IT=1,NT
                AT=1D0
                IF(IT.eq.1.or.IT.eq.NT) AT=5D-1
                SUMI1=SUMI1+AT*SOUMI0(IT,L)*BGMI1(1,IT,K)
  307 CONTINUE
              SUMI1=SUMI1*EPTA
      IF(IGEOM.ne.2) SUIPL(L)=SUIPL(L)+SUPL1
                     SUIMI(L)=SUIMI(L)+SUMI1
  303 CONTINUE
  301 CONTINUE
c internal angle averaged radiation field 
        DO 311 IX=1,II
              DO 313 ITAU=1,NT
              SUPL1=0D0
            IF(IGROUND.eq.2) THEN
                DO 315 IT=1,NT
                AT=1D0
                IF(IT.eq.1.or.IT.eq.NT) AT=5D-1
                SUPL1=SUPL1+AT*SOISTA(IT,IX)*BGMINA(IT,ITAU,IX)
  315 CONTINUE
            SUPL1=SUPL1*FILF1
            ENDIF
              IF(ITAU.ne.1) THEN
                DO 317 IT=1,ITAU
                AT=1D0
                IF(IT.eq.1.or.IT.eq.ITAU) AT=5D-1
                SUPL1=SUPL1+AT*SOISTA(IT,IX)*BGPLA(IT,ITAU,IX)
  317 CONTINUE
              ENDIF
              SUMI1=0D0
              IF(ITAU.ne.NT) THEN 
                DO 319 IT=ITAU,NT
                AT=1D0
                IF(IT.eq.ITAU.or.IT.eq.NT) AT=5D-1
                SUMI1=SUMI1+AT*SOISTA(IT,IX)*BGMIA(ITAU,IT,IX)
  319 CONTINUE
              ENDIF 
              DIMETA(IX,ITAU)=(SUPL1+SUMI1)*DELTAU/2D0 
  313 CONTINUE
  311 CONTINUE
      ENDIF 

      ISC=ISC+1 
      ENDDO 

c emergent radiation field through sides 
      IF(IGEOM.ne.0) THEN 
        DO 380 IX=1,II 
        X=XX(IX) 
        DO 380 IANG=1,NC
        K=IANG+(IX-1)*NC 
        L=K 
          ETA=UANG(IANG)
          EPTA=DELTAU/ETA 
          DNTETA=DNT*ETA 
          FLPL=0D0 
          FLMI=0D0 
          ITAU=1
          LTAU=NT
              DO 360 ITAU=1,NT
              LTAU=NT-ITAU+1
              ATAU=1D0
              IF(ITAU.eq.NT) ATAU=5D-1
         IF(IGEOM.eq.1.or.IGEOM.eq.2) THEN 
              SUPL2=0D0 
            IF(IGROUND.eq.2) THEN
                DO 362 IT=1,NT
                AT=1D0
                IF(IT.eq.1.or.IT.eq.NT) AT=5D-1
                SUPL2=SUPL2+AT*SOUMI(IT,L)*BGMEX(IT,ITAU,K)
  362 CONTINUE
            SUPL2=SUPL2*FILF1 
            ENDIF
                DO 364 IT=1,ITAU 
                AT=1D0 
                IF(IT.eq.1.or.IT.eq.ITAU) AT=5D-1 
                SUPL2=SUPL2+AT*SOUPL(IT,L)*BGPL2(IT,ITAU,K)
  364 CONTINUE 
              SUPL2=SUPL2*EPTA
              SUMI2=0D0 
                DO 366 IT=LTAU,NT 
                AT=1D0 
                IF(IT.eq.LTAU.or.IT.eq.NT) AT=5D-1 
                SUMI2=SUMI2+AT*SOUMI(IT,L)*BGMI2(LTAU,IT,K)
  366 CONTINUE 
              SUMI2=SUMI2*EPTA
c emergent side 
             FLPL=FLPL+ATAU*SUPL2 
             FLMI=FLMI+ATAU*SUMI2 
         ENDIF 
  360 CONTINUE
            FLPL=FLPL/DNTETA 
            FLMI=FLMI/DNTETA 
            SUIMIOUT(L)=SUIMIOUT(L)+FLMI 
              SUIPL(L)=SUIPL(L)+FLPL
              SUIMI(L)=SUIMI(L)+FLMI
  380 CONTINUE
      ENDIF 

      TAUOLD=TAU 
      YPHOLD=YPH 
      YOLD=Y
      FIL1OLD=FILF1 
      IGEOLD=IGEOM 
      RETURN
      END                  
c***********************************************************************


c***********************************************************************
C  calculation of  frequencies 
      SUBROUTINE FREQPOIN(XLMIN,XLMAX)
      IMPLICIT REAL*8(A-H,O-Z)   
      PARAMETER(IDEVC=12) 
      PARAMETER(MAXTAU=15)
      PARAMETER(MAXFRE=71,MAXANG=3,MAXSC=50)
      PARAMETER(II=MAXFRE,NSP=2*II,NC=MAXANG,ND=1,KK=II*NC,LL=ND*KK,
     1NLYR=MAXTAU,NL1=NLYR+1)
      PARAMETER(NEX=II+II*(II-1)/2,NCX=NC+NC*(NC-1)/2)
      PARAMETER(NX1=II+1)
      COMMON/JPWFRE/XX(II),XLOG(II),XKEV(II),EXXY(NEX),S0(II)
      COMMON/JPKN/TPH,TE,THOM,HX,HPHK02,NL,NT,NSCAT,NPOL 
      COMMON/QQX12/XX12(NX1) 
      DATA CKEV/0.0019569472D0/,EPS/1D-14/
      DATA NOUT/6/
      DLN=DLOG(1D1) 
          XL=XLMIN  
          XL1=XL-HX*0.5D0
          XX12(1)=DEXP(XL1*DLN)  
           DO 15 IX=1,II
           XLOG(IX)=XL 
           X=DEXP(XL*DLN)
           XL1=XL1+HX   
           XX12(IX+1)=DEXP(XL1*DLN) 
           XX(IX)=X                
           XEV=X/CKEV             
           XKEV(IX)=XEV          
           XL=XL+HX             
   15      CONTINUE            
      RETURN
      END  
c***********************************************************************
c***********************************************************************
C  calculation of cross-section S0(x) and axilliary quantities
      SUBROUTINE CROSCOMP
      IMPLICIT REAL*8(A-H,O-Z) 
      PARAMETER(IDEVC=12) 
      PARAMETER(MAXTAU=15)
      PARAMETER(MAXFRE=71,MAXANG=3,MAXSC=50)
      PARAMETER(II=MAXFRE,NSP=2*II,NC=MAXANG,ND=1,KK=II*NC,LL=ND*KK,
     1NLYR=MAXTAU,NL1=NLYR+1)
      PARAMETER(NEX=II+II*(II-1)/2,NCX=NC+NC*(NC-1)/2)
      PARAMETER(NX1=II+1)
      COMMON/JPWFRE/XX(II),XLOG(II),XKEV(II),EXXY(NEX),S0(II)
      COMMON/QWITCH/IGROUND,ICHAND,KEYRAY,IREDF,IDISTF,IGEOM,IPOL,ISP  
      COMMON/JPKN/TPH,TE,THOM,HX,HPHK02,NL,NT,NSCAT,NPOL 
      COMMON/QQX12/XX12(NX1) 
      COMMON /QQWYY/ YE,Y1,YR1,YR2,DY,YKY,DK02,DK2Y,G0 
      DATA CKEV/0.0019569472D0/,EPS/1D-14/
      DATA NOUT/6/
      DLN=DLOG(1D1)    
           DO 15 IX=1,II
           XL=XLOG(IX)
           X=XX(IX) 
            IF(IDISTF.eq.0) THEN 
c Maxwellian electrons
           CALL MXWELL(X,AB,XM,QM,PM,EPS) 
                DO 30 JX=1,IX 
                X1=XX(JX)
                ISX=JX+IX*(IX-1)/2 
                YXX=YE*(X-X1) 
                EXXY(ISX)=0D0
                IF(YXX.LT.7D1) EXXY(ISX)=DEXP(-YXX) 
   30           CONTINUE
            ELSEIF(IDISTF.eq.2) THEN
c Power-law electrons
           CALL PLCROS(X,AB,EPS) 
            ELSEIF(IDISTF.eq.4) THEN
c Maxwellian electrons + Power-law tail
           CALL PMCROS(X,AB,EPS) 
            ENDIF 
           S0(IX)=AB
   15      CONTINUE    
      RETURN
      END   
c***********************************************************************

c***********************************************************************
C angular distribution of soft radiation  
      REAL FUNCTION QANGLE*8(CI,ICHAND) 
      IMPLICIT REAL*8(A-H,O-Z)
      IF(ICHAND.EQ.0) QANGLE=1D0 
      IF(ICHAND.EQ.1) QANGLE=1D0+2.06D0*CI
      IF(ICHAND.EQ.2) QANGLE=1D0/CI
      RETURN
      END 
c***********************************************************************

c***********************************************************************
C calculation of cross-section for Power-law electron distribution 
c eq. (3.4.1) in Nagirner & Poutanen 1994 
      SUBROUTINE PLCROS(X,AB,EPS)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(NG=5) 
      COMMON /JPWZ/ Z,Z2,G,G2,GT,G1
      COMMON /JPPOW/PLIND,GMIN,GMAX,DNORPOW,DNORMAX,RATIO 
      COMMON/QQWG/UG(NG),AG(NG) 
      DATA PI/3.141592653589792384D0/ 
      DL21=DLOG(GMAX/GMIN)
c integral over electron distribution 
c loop over gamma 
      SUM=0D0 
      DO 10 I=1,NG
      DLU=UG(I)*DL21
      DLG=DLU+DLOG(GMIN)
      G=GMIN*DEXP(DLU)
      CALL COMWG
      GP=DEXP(-PLIND*DLG)
      GPZ=GP/Z*AG(I) 
c function of x and gamma (3.3.3) 
      SUM=SUM+PS0JP(X,EPS)*GPZ  
   10 CONTINUE 
      AB=SUM*DNORPOW*DL21*4D0*PI  
      RETURN
      END
c***********************************************************************
C calculation of cross-section for Maxwellian electrons + Power-law tail
c eq. (3.4.1) in Nagirner & Poutanen 1994 
      SUBROUTINE PMCROS(X,AB,EPS)
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(NG=5) 
      COMMON /JPWZ/ Z,Z2,G,G2,GT,G1
      COMMON /JPPOW/PLIND,GMIN,GMAX,DNORPOW,DNORMAX,RATIO 
      COMMON/QQWG/UG(NG),AG(NG) 
      COMMON /QQWYY/ YE,Y1,YR1,YR2,DY,YKY,DK02,DK2Y,G0 
      DATA PI/3.141592653589792384D0/ 
      DL21=DLOG(GMAX/GMIN)
      DLGM=DLOG(GMIN) 
c integral over power-law part of the electron distribution 
      SUM1=0D0 
      DO 10 I=1,NG
      DLU=UG(I)*DL21
      DLG=DLU+DLGM
      G=GMIN*DEXP(DLU)
      CALL COMWG
      GP=DEXP(-PLIND*DLG)
      GPZ=GP/Z*AG(I) 
c function of x and gamma (3.3.3) 
      SUM1=SUM1+PS0JP(X,EPS)*GPZ  
   10 CONTINUE 
      SUM1=SUM1*DNORPOW*DL21*4D0*PI  
c integral over Maxwellian part
            S=0D0
               DO 20 I=1,NG
               DGZ=UG(I)*DLGM
               G=DEXP(DGZ)
               DGY=(G-1D0)*YE
               EGA=0D0
               IF(DGY.lt.7d1) EGA=DEXP(-DGY)
               CALL COMWG
               S=S+PS0JP(X,EPS)*AG(I)*EGA*G
   20          CONTINUE
      S=S*DNORMAX*DLGM*YKY*YE
      AB=SUM1+S
      RETURN
      END
c***********************************************************************
C calculation of cross-section, mean energy of scattered photon, 
C dispersion of energy of scattered photon, and radiative pressure for   
C Maxwellian electron distribution. 
C THIS VERSION CALCULATES ONLY THE TOTAL  
C C S CROSS-SECTION:   AB (in units of Thomson 
C                                         cross-section);  
C INPUT: 
C 1. X - PHOTON ENERGY IN UNITS OF ELECTRON REST MASS m_ec^2 ; 
C 2. E - RELATIVE ERROR; 
C OUTPUT: 
C 3. AB - TOTAL CS CROSS-SECTION (in units of Thomson 
C         cross-section;  
C 4. XM - MEAN ENERGY OF SCATTERED PHOTON; 
C 5. QM - DISPERSION OF ENERGY OF SCATTERED PHOTON; 
C 6. PM - RADIATIVE PRESSURE;    
C  
      SUBROUTINE MXWELL(X,AB,XM,QM,PM,E)                              
      IMPLICIT REAL*8(A-H,O-Z)                                          
      COMMON /WAGUE/ UL(40),AL(40),AU(40),AE(40),N0 
      COMMON /QQWYY/ Y,YR,YR1,YR2,DY,YKY,DK02,DK2Y,G0 
      DATA C13 /0.3333333333333333333D0/,DE1/3.67879441171442D-01/       
c      X2=X*X                                   
      AB=0D0                                   
C      XM=0D0                             
C      QM=0D0                                   
C      PM=0D0                                   
      DO 600 L=1,N0                            
         U=UL(L)+1D0                           
         A=AU(L)                               
         AEL=AE(L)                             
         AGL=AL(L)                             
         G1=1D0+A*A*YR                         
         GY1=G1+YR                             
         YG1=G1*(G1+YR1)+YR2                   
         DZ1=DSQRT(1D0+G1)                     
         DR1=1D0/DZ1                           
         Z1=A*DZ1*DY                           
         GP1=G1+Z1                            
         ZP1=GP1*GP1                          
         GM1=1D0/GP1                          
         ZM1=GM1*GM1                          
c         CALL QSN(X*GP1,0,SSP,SSP1,SSP2,SP1,SP2,SP3,SP4,SP5,SP6,SP7,E) 
c         CALL QSN(X*GM1,0,SSM,SSM1,SSM2,SM1,SM2,SM3,SM4,SM5,SM6,SM7,E) 
         CALL QSN(X*GP1,0,SSP,E) 
         CALL QSN(X*GM1,0,SSM,E) 
         AB1=(ZP1*SSP+ZM1*SSM)*DR1                                    
C         XM1=((GP1*ZP1*(GY1*SP1+X*SP2)+GM1*ZM1*(GY1*SM1+X*SM2)))*DR1  
C         QM1=(ZP1*(0.5D0*SP6-GP1*(GY1*SP5+GP1*(SP7-YG1*SP4)))+        
C     +        ZM1*(0.5D0*SM6-GM1*(GY1*SM5+GM1*(SM7-YG1*SM4))))*DR1    
C         PM1=(ZP1*ZP1*SP1+ZM1*ZM1*SM1)*DR1                       
         G2=1D0+U*YR                                             
         GY2=G2+YR                                               
         YG2=G2*(G2+YR1)+YR2                                     
         DZ2=DSQRT(U*(1D0+G2))                                   
         DR2=1D0/DZ2                                             
         Z2=DZ2*DY                                               
         GP2=G2+Z2                                               
         ZP2=GP2*GP2                                             
         GM2=1D0/GP2                                             
         ZM2=GM2*GM2                                             
c         CALL QSN(X*GP2,0,SSP,SSP1,SSP2,SP1,SP2,SP3,SP4,SP5,SP6,SP7,E)
c         CALL QSN(X*GM2,0,SSM,SSM1,SSM2,SM1,SM2,SM3,SM4,SM5,SM6,SM7,E)
         CALL QSN(X*GP2,0,SSP,E) 
         CALL QSN(X*GM2,0,SSM,E) 
         AB2=(ZP2*SSP+ZM2*SSM)*DR2                                   
C         XM2=(GP2*ZP2*(GY2*SP1+X*SP2)+GM2*ZM2*(GY2*SM1+X*SM2))*DR2   
C         QM2=(ZP2*(0.5D0*SP6-GP2*(GY2*SP5+GP2*(SP7-YG2*SP4)))+       
C     +        ZM2*(0.5D0*SM6-GM2*(GY2*SM5+GM2*(SM7-YG2*SM4))))*DR2   
C         PM2=(ZP2*ZP2*SP1+ZM2*ZM2*SM1)*DR2 
         AB=AB+AGL*(AEL*AB1+DE1*AB2) 
C         XM=XM+AGL*(AEL*XM1+DE1*XM2) 
C         QM=QM+AGL*(AEL*QM1+DE1*QM2) 
C         PM=PM+AGL*(AEL*PM1+DE1*PM2) 
  600 CONTINUE  
      AB=AB*DK2Y                     
C      XM=XM*DK2Y                     
C      QM=QM*DK2Y                     
C      PM=PM*DK2Y                     
      RETURN                         
      END                            


C  S- functions. if K=0, then only  SS -- mean cross-section; 
C if  K=1, then  SS1 &  S1 -- mean frequency;
C if K=2, all functions -- and the square of frequency
C x -- small
c      SUBROUTINE QSN(X,K,SS,SS1,SS2,S1,S2,S3,S4,S5,S6,S7,E) 
      SUBROUTINE QSN(X,K,SS,E) 
      IMPLICIT REAL*8(A-H,O-Z)                             
      DATA XT /4D-1/                                       
      DATA C13/0.333333333333333D0/                        
      IF (X.LT.XT) GO TO 110                               
      X0=1D0/(1D0+2D0*X)                                   
      X2=X*X                                               
      XR=1D0/X                                             
      DLX0=-DLOG(X0)                                       
      X02=X0*X0                                            
      X2R=XR*XR                                            
      SS=0.375D0*(4D0-(2D0*XR+2D0-X)*DLX0+2D0*X2*(1D0+X)*X02)*X2R    
      IF(K.EQ.0) RETURN                                              
      SS1=0.125D0*X2R*XR*(3D0*DLX0+4D0*X2-4.5D0*X-X*X0*(1.5D0+X*X02))
      S1=(SS-SS1)*XR                             
      S2=(SS1-S1)*XR                             
      IF(K.EQ.1) RETURN                          
            SS2=(1D0+X*(4D0+X*(7D0+4.5D0*X)))*X02*X02
      S3=(SS1-SS2)*XR                                
            S4=(S1-S3)*XR                            
      GO TO 160                                      
  110 A=1D0                                          
      B=-2D0*X                                       
      DN=0D0                                         
      SS=0D0                                         
      IF (K.EQ.0) GO TO 120                          
      SS1=0D0                                        
      S1=0D0                                         
      S2=0D0                                         
      IF (K.EQ.1) GO TO 120                          
      SS2=0D0                                        
      S3=0D0                                         
      S4=0D0                                         
  120 DN1=DN+1D0                                     
      DN2=DN+2D0                                     
      DN3=DN+3D0                                     
      D1=1D0/DN1                                     
      D2=1D0/DN2                                     
      D3=1D0/DN3                                     
      DN4=DN+4D0                                     
      DN23=DN2*DN3                                   
      DN234=0.25D0*DN23*DN4                          
      SS=SS+A*(DN2+2D0*D1+8D0*D2-16D0*D3)            
      IF (K.EQ.0) GO TO 140                          
            D4=1D0/DN4                               
      D5=1D0/(DN4+1D0)                               
            DN13=DN1*DN3                             
                  DN14=DN1*DN4                       
      SS1=SS1+A*(DN23-6D0+24D0*D3)                   
      S1=S1+A*(DN13-6D0-6D0*D2-24D0*D3+72D0*D4)      
      S2=S2+A*(DN14-6D0-12D0*D3-72D0*D4+144D0*D5)    
      IF (K.EQ.1) GO TO 140                          
      DN134=0.25D0*DN13*DN4                          
      DN124=0.25D0*DN2*DN14                          
      SS2=SS2+A*(DN234+2D0-DN)                       
      S3=S3+A*(DN134+7D0-DN-24D0*D4)                 
      S4=S4+A*(DN124+12D0-DN+6D0*D3+24D0*D4-96D0*D5) 
  140 IF (DABS(4D0*A)*DN234.LT.E) GO TO 150          
      DN=DN+1D0                                      
      A=A*B                                          
      GO TO 120                                      
  150 SS=SS*0.375D0                                  
      IF (K.EQ.0) RETURN                             
      SS1=SS1*0.125D0                                
      S1=S1*0.25D0                                   
      S2=S2*0.25D0                                   
      IF (K.EQ.1) RETURN                             
            SS2=SS2*0.125D0                          
      S3=S3*0.25D0                                   
      S4=S4*0.5D0                                    
  160 S5=3D0*S4-4D0*S3                               
      S7=S3-0.5D0*S4                                 
      S6=2D0*SS2-6D0*S7                              
      RETURN                                         
      END                                            


      SUBROUTINE COMWY(EPS)                                           
C*** CALCULATES SOME FUNCTION OF ELECTRON TEMPERATURE. 
C* * INPUT: 
C* * 1. Y - INVERSE TEMPERATURE (m_c^2/kT_e) .
C* * OUTPUT: 
C* * 2. YR= 1/Y, YR1=2/Y, YR2=2/(Y*Y), DY=1./SQRT(Y)  
C* * 3. YKY=EXP(-Y)/K_2(Y) 
C* * 4. DK2Y=YKY/(2.*SQRT(Y))
C* * 5. DK02=K_0(Y)/K_2(Y) 
C* * 6. DK12=K_1(Y)/K_2(Y) 
C* * 7. GAMMEAN=3D0/Y + K_1(y)/K_2(y) 
C* * 8. GAMSQUA=1D0+3/y*[K_3(y)/K_2(y)]  
C***
      IMPLICIT REAL*8 (A-H,O-Z)                                       
      PARAMETER(NG=5) 
      COMMON /QQWYY/ Y,YR,YR1,YR2,DY,YKY,DK02,DK2Y,G0 
      COMMON /JPWYAV/ GAMMEAN,ZMEAN,GAMSQUA,DK12 
      COMMON /JPPOW/PLIND,GMIN,GMAX,DNORPOW,DNORMAX,RATIO 
      COMMON /JPWZ/ Z,Z2,G,G2,GT,G1  
      COMMON/QQWG/UG(NG),AG(NG)
      COMMON/QWITCH/IGROUND,ICHAND,KEYRAY,IREDF,IDISTF,IGEOM,IPOL,ISP  
      DIMENSION DK(20)                                              
      DATA PI/3.141592653589792384D0/,C13/0.33333333333333333D0/
C    PIQR=1./SQRT(2.*PI)
      DATA PIQR/0.39894228040143267794D0/
c RATIO - ratio of the energy in power-law tail to the Maxwellian 
c         part of the electron distribution 
c CUT*<g> - energy of Maxwellian cut-off
c      DATA CUT/2D0/  
      DATA DIFG/1D-5/  
      IF(IDISTF.eq.0.or.IDISTF.eq.3.or.IDISTF.eq.4) THEN 
      PI2=2D0/PI                                                    
      CALL QDBESK(Y,4,DK,8D0,6)                                      
      DK2=DK(3) 
      YR=1D0/Y                
      G0=-YR*DLOG(EPS)        
      IF(Y.LT.1D0) G0=-YR*DLOG(EPS*G0) 
      DK02=DK(1)/DK2
      DK12=DK(2)/DK2
      GAMMEAN=3D0*YR + DK12 
      ZMEAN=DSQRT(GAMMEAN*GAMMEAN-1D0)
      GAM1=GAMMEAN-1D0 
      GAMSQUA=1D0+3D0/Y*DK(4)/DK2 
      YR1=2D0*YR             
      YR2=YR1*YR             
      DY=DSQRT(YR)           
      DK2R=1D0/DK2            
        IF(Y.GT.8D0) THEN 
        YPI=DSQRT(PI2*Y)        
        YKY=YPI/DK2           
        DK2Y=PIQR*DK2R          
        ELSE 
        EDY=DEXP(-Y)            
        YKY=EDY/DK2      
        DK2Y=0.5D0*EDY*DY*DK2R  
        ENDIF 
      ENDIF 

       IF(IDISTF.eq.2.or.IDISTF.eq.4.and.GMIN.gt.GMAX) THEN 
       GMID=GMIN
       GMIN=DMIN1(GMIN,GMAX)
       GMAX=DMAX1(GMID,GMAX)
       ENDIF 
       IF(IDISTF.eq.4.and.GMIN.eq.GMAX) IDISTF=3 

       IF(IDISTF.eq.3.or.IDISTF.eq.4) THEN 
c       GMIN=CUT*GAMMEAN 
       DLGM=DLOG(GMIN)
            SUM1=0D0
            SUM3=0D0
               DO 10 I=1,NG
               DGZ=UG(I)*DLGM
               G=DEXP(DGZ)
               DGY=(G-1D0)*Y
               EGA=0D0
               IF(DGY.lt.7d1) EGA=DEXP(-DGY)
               CALL COMWG
               SUM1=SUM1+Z*G2*AG(I)*EGA
       IF(IDISTF.eq.4) SUM3=SUM3+Z*G2*G*AG(I)*EGA
   10  CONTINUE
       SUM1=SUM1*DLGM*Y*YKY
       SUM3=SUM3*DLGM*Y*YKY
       ENDIF

c normalization constant of power-law part 
       IF(IDISTF.eq.2.or.IDISTF.eq.4) THEN 
        DLGMAMI=DLOG(GMAX/GMIN)
        P1=PLIND-1D0
        PI4=4D0*PI
        DNORPOW=1D0/(PI4*DLGMAMI)
        IF(DABS(P1).GT.1D-5) DNORPOW=P1/PI4*DEXP(P1*DLOG(GMIN))/
     *  (1D0-DEXP(-P1*DLGMAMI))
       ENDIF
 
c normalization constant for maxwellian part  
       IF(IDISTF.eq.3) THEN 
       DNORMAX=1D0/SUM1 
       ENDIF

       IF(IDISTF.eq.4) THEN 
       DLGMAMI=DLOG(GMAX/GMIN)
       GMIN2=GMIN*GMIN
       GMIN3=GMIN*GMIN2 
       ZMIN=DSQRT(GMIN*GMIN-1D0) 
       GMY=Y*(GMIN-1D0) 
       EGA=0D0 
       IF(GMY.lt.7d1) EGA=DEXP(-GMY) 
       YKYEGA=YKY*EGA
       DGPMI=(1D0+PLIND)*DLOG(GMIN) 
       EGPMI=0D0
       IF(DGPMI.lt.5D1) EGPMI=DEXP(DGPMI)
       DNORPOW=Y*YKYEGA*ZMIN*EGPMI/(PI*4D0)/DNORPOW
       DNORMAX=1D0/(SUM1+DNORPOW)
       DNORPOW=DNORMAX*DNORPOW 
       SUM4=PI4*DLGMAMI
       P2=PLIND-2D0
       IF(DABS(P2).GT.1D-5) SUM4=PI4*DEXP(-P2*DLOG(GMIN))*
     * (1D0-DEXP(-P2*DLGMAMI))/P2
       RATIO=DNORPOW*SUM4/(DNORMAX*SUM3) 
       ENDIF

      RETURN   
      END      


C
C AUXILARY QUANTITIES DEPENDING ON X,X1 AND MU ; 
C NEEDED TO CALL FIRST  COMX(X,X1) SUBROUTINE 
C 
      SUBROUTINE COMXM(QM)                                 
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON /JPWX/ XX,XX1,SX,DX,T,TR,XM2,DX2,GM1,GA1,RB,RWB,NR
      COMMON /JPWXM/ GST,GST1,QO,QR,QOT,Q2,Q,QT,QMT,PM,RS,QS
      DATA EPS/1D-10/ 
      QO=T*QM    
      Q2=DX2+2D0*QO                   
      Q=DSQRT(Q2)                     
      IF(QM.gt.EPS) THEN 
      QR=1D0/QO                       
      QOT=2D0*QR                      
      QT=2D0/Q                        
      QMT=2D0/QM                      
      PM=2D0-QM                       
      RS=PM/QM                        
      QS=T*PM                         
      GST=(DX+Q*DSQRT(1D0+QOT))*5D-1
      GST1=GST-1D0
      ELSE 
      GST=1D0 
      GST1=0D0
      QT=0D0 
      ENDIF 
      RETURN                                          
      END                           


C 
C AUXILARY QUANTITIES DEPENDING ON GAMMA 
C 
      SUBROUTINE COMWG            
      IMPLICIT REAL*8(A-H,O-Z)     
      COMMON /JPWZ/ Z,Z2,G,G2,GT,G1  
      G1=G-1D0                     
      G2=G*G                       
      GT=2D0*G                     
      Z2=G1*(G1+2D0)               
      Z=DSQRT(Z2)                  
      RETURN                       
      END                          

      SUBROUTINE COMWZ
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON /JPWZ/ Z,Z2,G,G2,GT,G1
      Z2=Z*Z
      G2=Z2+1D0
      G=DSQRT(G2)
      G1=G-1D0
      GT=2D0*G
      RETURN
      END

C 
C AUXILARY QUANTITIES DEPENDING ON  X AND X1
C 
      SUBROUTINE COMX(X,X1)                                 
      IMPLICIT REAL*8(A-H,O-Z)                              
      COMMON /JPWX/XX,XX1,SX,DX,T,TR,XM2,DX2,GM1,GA1,RB,RWB,NR
      COMMON /JPWFT/TRT,TRF,TRR,TRS,TR2,XS,XS1,XA,XA1,XXA,XXA1
      COMMON /JPWFW/TRO,TT,TOO,TTT                            
      DATA C13/0.33333333333333D0/,C23/0.66666666666666D0/  
      XS=1D0+2D0*X                                          
      XA=X/XS                                               
      XS1=1D0+2D0*X1                                        
      XA1=X1/XS1                                            
      XX=2D0*XA*X     
      XX1=-2D0*XA1*X1 
      SX=X+X1         
      DX=X-X1         
      T=X*X1          
      TT=T*T          
      TR=1D0/T        
      XM2=SX*DX       
      DX2=DX*DX       
      TR2=TR*TR       
      TRT=-2D0*TR     
      TRF=4D0*TR      
      TRR=TR+1D0      
      TRS=5D-1*TR2                    
      TTT=DX2*TR2-TRF                 
      TOO=1D0-TRT                     
      XMAX=X                          
      IF(X.LT.X1) XMAX=X1             
      XMIN=X                          
      IF(X.GT.X1) XMIN=X1             
      RB=4D0/XMAX                     
      RWB=RB-4D0*C13*XMIN/(XMAX*XMAX) 
      GA=0.5D0*(DX+SX*DSQRT(TRR))     
      XXA=X-XA1                       
      XXA1=X1-XA                      
      XXX=(1D0-2D0*X)/XS              
      GA1=GA-1D0                      
      IF(XXA.GT.-0.1D0)      GA1=0.25D0*XS1*XS1*TR*XXA*XXA/(GA+1D0-DX)  
      IF(XXA1.GT.-0.1D0*XXX) GA1=0.25D0*XS*XS*TR*XXA1*XXA1/(GA+1D0)+DX  
      IF(XXA.LE.0D0) NR=1                 
      IF(XXA1.LE.0D0) NR=2                
      IF((XXA.GT.0D0).AND.(X.LE.X1)) NR=3 
      IF((XXA1.GT.0D0).AND.(X1.LT.X)) NR=4
      IF(NR-3) 10,20,30                   
   10 GM1=GA1                             
      GO TO 40                            
   20 CONTINUE                            
      GM1=0D0                             
      GO TO 40                           
   30 CONTINUE                             
      GM1=DX                             
   40 CONTINUE                           
      RETURN                             
      END                                


C 
C REDISTRIBUTION MATRIX : S(X,X1,MU)
C INPUT: 
C 1. X, X1 - ENERGIES OF SCATTERED AND  INITIAL PHOTONS, RESPECTIVELY; 
C 2. QM=1-MU , AND MU - THE COSINE OF THE SCATTERING ANGLE; 
C OUTPUT: 
C 3. S, SI, SQ, SU, SV, SA=SQ-SU  (X,X1,MU) 
C 
      SUBROUTINE REDFUNC(X,X1,QM,S) 
      IMPLICIT REAL*8(A-H,O-Z)                                    
      COMMON /JPWXM/ GST,GST1,QO,QR,QOT,Q2,Q,QT,QMT,PM,RS,QS 
      COMMON /JPWX/ XX,XX1,SX,DX,T,TR,XM2,DX2,GM1,GA1,RB,RWB,NR
      COMMON /QQWYY/ YE,Y1,YR1,YR2,DY,YKY,DK02,DK2Y,G0 
      COMMON /JPWYAV/ GAMMEAN,ZMEAN,GAMSQUA,DK12 
      COMMON /JPPOW/PLIND,GMIN,GMAX,DNORPOW,DNORMAX,RATIO 
      COMMON/QWITCH/IGROUND,ICHAND,KEYRAY,IREDF,IDISTF,IGEOM,IPOL,ISP  
      COMMON /JPWZ/ Z,Z2,G,G2,GT,G1   
C PI332=3/(32*PI)      
      DATA PI/3.141592653589792384D0/,C13/0.33333333333333333D0/
      DATA PI332/0.2984155182973039D-01/,EPS/1D-10/ 
      DATA NOUT/6/
      CALL COMXM(QM) 
c Maxwellian electron distribution 
          IF(IDISTF.eq.0.or.IDISTF.eq.3.or.IDISTF.eq.4) THEN 
            EXGST=0D0                                            
            YGST=YE*GST1                                    
            IF(YGST.LT.7D1) EXGST=DEXP(-YGST)     
            CONST=PI332*YKY*EXGST
          ENDIF 
c exact Compton redistribution function for Maxwellian electrons 
        IF(IDISTF.eq.0) THEN 
          IF(IREDF.eq.0) THEN     
            CALL MAXGALA(X,X1,QM,S) 
            S=S*CONST    
          ENDIF 
c isotropic approximation in electron rest frame 
c S=1/8*pi exp(-y g_*)/K_2(y)  g_* exact=(x-x_1+Q*sqrt(1+2/q))/2
          IF(IREDF.eq.1) THEN 
          S=CONST*C13*2D0*QT   
          ENDIF 
        ENDIF
c Maxwellian electron distribution with cut-off 
c and Maxwellian part of RF for Maxw + power-law
          IF(IDISTF.eq.3.or.IDISTF.eq.4) THEN 
           IF(GMIN.gt.GST) THEN 
c exact RF 
            IF(IREDF.eq.0) CALL MAXCUT(X,X1,QM,SM)  
c isotropic approximation 
            IF(IREDF.eq.1) THEN 
            DG=YE*(GMIN-GST) 
            EGA=0D0 
            IF(DG.lt.7d1) EGA=DEXP(-DG) 
            SM=C13*2D0*QT*(1D0-EGA)*DNORMAX  
            ENDIF
            SM=SM*CONST 
           ELSE 
            SM=0D0 
           ENDIF
          ENDIF
          IF(IDISTF.eq.3) S=SM 
c power-law part of 
c Maxwellian electron distribution with power-law after CUT*GAMMEAN 
          IF(IDISTF.eq.4) THEN 
           IF(GMAX.gt.GST) THEN 
c exact RF 
            IF(IREDF.eq.0) CALL POWMATR(X,X1,QM,SP)
c isotropic approximation using fact that CUT*GAMMEAN>>1 
            IF(IREDF.eq.1) THEN 
            GM=DMAX1(GMIN,GST) 
            DLGMM=DLOG(GMAX/GM) 
            P1=PLIND+1D0 
            P3=PLIND+3D0 
            GP3=DEXP(-P3*DLOG(GM)) 
            GMM1=DEXP(-P1*DLGMM) 
            GMM3=DEXP(-P3*DLGMM) 
            SP=0.25D0*QT*DNORPOW*GP3*(GM*GM*(1D0-GMM1)/P1+
     1      5D-1*(1D0-GMM3)/P3) 
            ENDIF
           S=SM+SP 
           ELSE 
           S=0D0 
           ENDIF 
          ENDIF                                          
c monoenergetic electrons 
         IF(IDISTF.eq.1) THEN
c exact 
          IF(IREDF.eq.0) THEN 
            IF(GAMMEAN .ge. GST) THEN 
            G=GAMMEAN 
            CALL COMWG
            CALL RDSTRB(X,X1,QM,REDF)
            S=PI332/(ZMEAN*GAMMEAN)*REDF
            ELSE
            S=0D0 
            ENDIF
          ENDIF 
c isotropic approximation in electron rest frame 
          IF(IREDF.eq.1)  THEN 
           IF(GAMMEAN .ge. GST) THEN  
           S=0.125D0*QT/(2D0*PI*GAMMEAN*ZMEAN)
           ELSE 
           S=0D0 
           ENDIF 
          ENDIF 
         ENDIF 
c power-law electron distribution 
c IF IREDF.eq.0 then exact
c IF IREDF.eq.1 then R=4/(3*Q) - isotropic approximation  
         IF(IDISTF.eq.2) CALL POWMATR(X,X1,QM,S) 
      RETURN                                               
      END                                                  


C 
C SCATTERING MATRIX FOR MAXWELLIAN ELECTRONS 
C 
      SUBROUTINE MAXGALA(X,X1,QM,S) 
      IMPLICIT REAL*8(A-H,O-Z)                                    
      COMMON /WAGUE/ UL(40),AL(40),AU(40),AE(40),N0        
      COMMON /QQWYY/ YE,Y1,YR1,YR2,DY,YKY,DK02,DK2Y,G0 
      COMMON /JPWXM/ GST,GST1,QO,QR,QOT,Q2,Q,QT,QMT,PM,RS,QS
      COMMON /JPWZ/ Z,Z2,G,G2,GT,G1   
            S=0D0
               DO 10 I=1,N0
               DGZ=UL(I)*Y1 
               G=DGZ+GST
               CALL COMWG
               CALL RDSTRB(X,X1,QM,R) 
               S=S+R*AL(I)                                       
   10          CONTINUE                                          
      RETURN                                               
      END                                                  

C SCATTERING MATRIX FOR MAXWELLIAN ELECTRONS with cut-off on 
C CUT*GAMMEAN energy 
C 
      SUBROUTINE MAXCUT(X,X1,QM,S)
      IMPLICIT REAL*8(A-H,O-Z)                                    
      PARAMETER(NG=5) 
      COMMON /JPPOW/PLIND,GMIN,GMAX,DNORPOW,DNORMAX,RATIO 
      COMMON/QWITCH/IGROUND,ICHAND,KEYRAY,IREDF,IDISTF,IGEOM,IPOL,ISP  
      COMMON/QQWG/UG(NG),AG(NG)
      COMMON /QQWYY/ YE,Y1,YR1,YR2,DY,YKY,DK02,DK2Y,G0 
      COMMON /JPWXM/ GST,GST1,QO,QR,QOT,Q2,Q,QT,QMT,PM,RS,QS
      COMMON /JPWZ/ Z,Z2,G,G2,GT,G1   
            S=0D0
        IF(GMIN.gt.GST) THEN 
        DGMS=DLOG(GMIN/GST) 
        DGST=DLOG(GST) 
               DO 10 I=1,NG 
               DGZ=UG(I)*DGMS+DGST 
               G=DEXP(DGZ) 
               DGY=(G-GST)*YE 
               EGA=0D0 
               IF(DGY.lt.7d1) EGA=DEXP(-DGY) 
               CALL COMWG 
               CALL RDSTRB(X,X1,QM,R) 
               S=S+R*AG(I)*EGA*G 
   10          CONTINUE 
      S=S*DNORMAX*DGMS*YE 
      ENDIF  
      RETURN                                               
      END 
                                                 

C Procedure calculates the S(x,x1,mu) matrix for the power-law electron 
C distribution:
C P- power-law index; GMIN,GMAX- limits of distribution; 
C X,X1 - frequencies in units of m_ec2; QM=1-mu; 
C GST=gamma_*; 
C need in advance calling the COMXX(X,X1) prosedure
      SUBROUTINE POWMATR(X,X1,QM,S)
      IMPLICIT REAL*8(A-H,O-Z)                 
      PARAMETER(NG=5) 
      COMMON /JPWZ/ Z,Z2,G,G2,GT,G1             
      COMMON/QQWG/UG(NG),AG(NG)
      COMMON /JPPOW/PLIND,GMIN,GMAX,DNORPOW,DNORMAX,RATIO 
      COMMON/QWITCH/IGROUND,ICHAND,KEYRAY,IREDF,IDISTF,IGEOM,IPOL,ISP  
      COMMON /JPWXM/ GST,GST1,QO,QR,QOT,Q2,Q,QT,QMT,PM,RS,QS
      DATA C38/0.375D0/,C13/0.3333333333333333D0/
      DATA C23/0.6666666666666667D0/
      DATA PI /3.141592653589792384D0/,GCR/2D0/ 
      S=1D-40
      IF(GST.GE.GMAX) RETURN
      IF(IREDF.eq.1) R=C23*QT
      GM=DMAX1(GST,GMIN)

       SUM1=0D0
       if(GCR.lt.GMAX) then
c int d(ln g)/z g^-p R(g)
      GMM=DMAX1(GM,GCR)
      DL21=DLOG(GMAX/GMM)
      DL1=DLOG(GMM)
      DO 10 I=1,NG
      AGI=AG(I)
      DLU=UG(I)*DL21
      DLG=UG(I)*DL21+DL1
      G=DEXP(DLG)
      CALL COMWG
      GP=DEXP(-PLIND*DLG)
      GPZ=GP/Z*AGI
      IF(IREDF.eq.0) CALL RDSTRB(X,X1,QM,R)
      SUM1=SUM1+R*GPZ
   10 CONTINUE
      SUM1=SUM1*DL21
       endif

       SUM2=0D0
       if(GCR.gt.GM) then
c int dz/g^2 g^-p R(g)
       GMM=DMIN1(GCR,GMAX)
       ZMM=DSQRT(GMM*GMM-1D0)
       ZM=DSQRT(GM*GM-1D0)
       DL21=ZMM-ZM
      DO 20 I=1,NG
      AGI=AG(I)
      Z=UG(I)*DL21+ZM
      CALL COMWZ
      DLG=DLOG(G)
      GP1=DEXP(-(PLIND+2D0)*DLG)
      IF(IREDF.eq.0) CALL RDSTRB(X,X1,QM,R)
      SUM2=SUM2+R*GP1*AGI
   20 CONTINUE
      SUM2=SUM2*DL21
       endif

      S=SUM1+SUM2
      S=S*C38*DNORPOW
      RETURN  
      END    
c ****end of POWMATR 

C 
C SCATTERING MATRIX R(X,X1,MU,GAMMA) 
C needed first to call subroutines: 
C COMX(X,X1); COMXM(QM); COMWG 
C  INPUT: 
C  1. X,X1 - ENERGIES OF  SCATTERED AND INITIAL PHOTONS;  
C  2. QM=1-MU (MU- COSINE OF SCATTERING ANGLE); 
C  3.  GAMMA=G (in COMMON/WZ/) 
C  OUTPUT: R,RI,RQ,RU,RV, RF=RQ-RU 
C 
      SUBROUTINE RDSTRB(X,X1,QM,R) 
      IMPLICIT REAL*8(A-H,O-Z)                               
      COMMON /JPWZ/ Z,Z2,G,G2,GT,G1                            
      COMMON /JPWX/ XX,XX1,SX,DX,T,TR,XM2,DX2,GM1,GA1,RB,RWB,NR
      COMMON /JPWXM/ GST,GST1,QO,QR,QOT,Q2,Q,QT,QMT,PM,RS,QS
      DATA PI /3.141592653589792384D0/,EPS/1D-10/                        
      IF(Z.LE.0D0)  stop 
      IF(QM.GT.2D0) stop 
      IF(X.LE.0D0)  stop
      IF(X1.LE.0D0) stop
      IF(QM.lt.EPS) THEN 
       R=0D0 
      RETURN  
      ELSE 
      A=DSQRT(Z2-(GT-X)*X+QMT)        
      A1=DSQRT(Z2+(GT+X1)*X1+QMT)     
      U=SX*(GT-DX)/(A+A1)             
      U2=U*U                                  
      V=A*A1                                  
      VR=1D0/V                                
      VQ=QR*VR                                
      UV=U*VR                                 
      UVQ=UV*QR                               
      QU=U+Q                                  
      UQ=U-Q                                   
      IF(UQ.LT.3D-1) GO TO 500               
      QSR=1D0/QS
      UQQ=U2-Q2                               
      RF=0.5D0*UV*VQ*VQ*(U2*UQQ+V*(5D0*U2-3D0*Q2)) 
      RG=QT+UV*(1D0-QOT)                      
      R=RF+RG                                 
      RETURN                                  
  500 GX=G-DX                                 
      DG=G-GST
      C=2D0/(G*GX+RS+T-QO+V)                  
      D=2D0*(GX+GST)*DG                        
      CD=C*D                                  
      UQQ=QS*CD                               
      QQQ=CD/QU                               
      RFO=0.5D0*UV*VQ*VQ*UQQ*(U2+5D0*V)       
      RGO=QT+UV                                
      RHO=UVQ*CD                              
      DFG=UV*VQ*((1D0-QS*C)*D+RS*(2D0+QO))    
      R=RFO+RGO-DFG                           
      ENDIF 
      RETURN                                 
      END                                                              


c***********************************************************************
c**************************************************
c function Psi_0 (x)*G*Z used in calculation of the total
c cross-section for Compton scattering
c eq. (3.3.3) in Nagirner & Poutanen 1994 
      FUNCTION PS0JP(X,E) 
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON /JPWZ/ Z,Z2,G,G2,GT,G1
      GZP=G+Z 
      GZP2=GZP**2  
      ARG1=X*GZP 
      ARG2=X/GZP 
      PS0JP=(PSIJP(ARG1,E)*GZP2-PSIJP(ARG2,E)/GZP2)/4D0  
      RETURN
      END
c**************************************************
c function psi_10 (x) used in calculation of the total 
c cross-section for Compton scattering 
c eq. (3.3.10) in Nagirner & Poutanen 1994 
      FUNCTION PSIJP(X,E) 
      IMPLICIT REAL*8(A-H,O-Z) 
      DATA XT /4D-1/
      DATA C13,C43 /0.3333333333333333D0,1.3333333333333333D0/  
      IF(X.LT.XT) GO TO 210 
      X0=1D0/(1D0+2D0*X)   
      X02=X0*X0           
      X03=X0*X02         
      X2=X*X            
      XR=1D0/X         
      DLX0=-DLOG(X0)  
      G0=GINT(X,E)      
      XR2=XR*XR     
      P10=0.75D0*XR2*((X+4.5D0+2D0*XR)*DLX0-4D0-X+X2*X0-2D0*G0)       
      GO TO 240 
  210 DN=0D0    
      A=1D0     
      B=-2D0*X  
      P10=0D0      
  220 DN1=DN+1D0 
      DN2=DN+2D0
      DN3=DN+3D0  
      DNQ=DN*DN    
      D1=1D0/DN1     
      D2=1D0/DN2    
      D3=1D0/DN3      
      D33=D3*D3         
      P10=P10+A*(1D0+D2*(2D0*D1+8D0*D2-16D0*D3)) 
      IF(DABS(A)*(DNQ+1D0).LT.E) GO TO 230  
      A=B*A                                
      DN=DN+1D0     
      GO TO 220    
  230 CONTINUE 
      P10=P10*0.75D0
  240 PSIJP=P10 
      RETURN  
      END  

      FUNCTION GINT(X,E)    
      IMPLICIT REAL*8(A-H,O-Z) 
      DATA XT,XTR /0.4D0,0.625D0/ 
      DATA PI,DL/3.141592653589792384D0,0.6931471805599453094D0/
      DATA PIH,PITW /1.64493406684823D0,0.822467033424113D0/   
      IF (X.LT.XT) GO TO 40                              
      IF (X.GT.XTR) GO TO 40  
      GINT=PITW                  
      IF (X.EQ.0.5D0) RETURN  
      Y=2D0*X                 
      Z=DLOG(Y)               
      GINT=GINT+DL*Z                
      IF (X.GT.0.5D0) GINT=GINT+0.5D0*Z*Z
      IF (X.GT.0.5D0) Y=1D0/Y  
      Y=1D0-Y   
      B=1D0    
      U=0D0     
      DK=1D0    
      V=0D0     
      A=1D0     
   30 A=0.5D0*A 
      V=V+A/DK  
      DK=DK+1D0 
      U=U+B*V/DK  
      B=B*Y       
      IF (B.GT.E) GO TO 30   
      U=U*Y*Y                
      IF (X.GT.0.5D0) GINT=GINT-U  
      IF (X.LT.0.5D0) GINT=GINT+U  
      RETURN                 
   40 IF (X.LT.XT) Y=2D0*X   
      IF (X.GT.XTR) Y=0.5D0/X
      GINT=1D0                 
      A=-Y                 
      DN=1D0     
      B=1D0      
   50 DN=DN+1D0  
      B=A*B      
      GINT=GINT+B/(DN*DN)  
      IF (DABS(B).GT.E) GO TO 50 
      GINT=Y*GINT          
      IF (X.LT.XT) RETURN  
      Z=DLOG(Y)     
      GINT=PIH+0.5D0*Z*Z-GINT     
      RETURN       
      END         
c*************************************************************************
c subroutine for computing photon supply from outside of the active region
c temperature distribution is assumed as follows: 
c T=Tin for r < Rc
c T=Tin * r**(-3/4) for r> Rc
      SUBROUTINE OUTFLUX(TAUVEC,DOU,TAU,T0TAUG) 
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(NG=5)
      PARAMETER(NCA=3)
      PARAMETER(MAXTAU=15)
      PARAMETER(MAXFRE=71,MAXANG=3,MAXSC=50)
      PARAMETER(II=MAXFRE,NSP=2*II,NC=MAXANG,ND=1,KK=II*NC,LL=ND*KK,
     1NLYR=MAXTAU,NL1=NLYR+1)
      PARAMETER(NEX=II+II*(II-1)/2,NCX=NC+NC*(NC-1)/2)
      PARAMETER(NX1=II+1)
      COMMON/JPWFRE/XX(II),XLOG(II),XKEV(II),EXXY(NEX),S0(II) 
      COMMON/QWITCH/IGROUND,ICHAND,KEYRAY,IREDF,IDISTF,IGEOM,IPOL,ISP
      COMMON/JPKN/TPH,TE,THOM,HX,HPHK02,NL,NT,NSCAT,NPOL
      COMMON/QQCM2/A(II),UANG(NC),AANG(NC),AC(LL),AINT(LL),CINT(LL) 
      DIMENSION TAUVEC(MAXTAU),DOU(MAXTAU,KK)
      DATA CKEV/0.0019569472D0/,PI/3.141592653589792384D0/,EPS/1D-14/
      DATA C13/0.33333333333333333D0/,C12/5D-1/,CMIN/7D1/
      TIN=TPH  

c beginning of tau loop
       sum=0d0
       do itau=1,nt
       sum1=0d0
c height in units of the radius
       taut=tauvec(itau)*t0taug 
       taut2=taut**2 
       if(IGEOM.eq.1) rr=1d0
       if(IGEOM.eq.2) then 
        if(dabs(taut2-1d0).gt.eps) then 
        rr=dsqrt(1d0-taut2) 
        else 
        rr=0d0
        endif 
       endif 
      
       dmu1=dsqrt((1d0+rr)/2d0)
       dmu2=taut/dsqrt(1d0+taut2) 
       dmu3=taut/dsqrt(taut2+(1d0+rr)**2) 
       
c beginning of mu loop
         do j=1,nc 
         dmu=uang(j)
         dms=dsqrt(1d0-dmu**2)
       t=taut*dms/dmu  
       t2=t*t 
c beginning of x loop         
       sum2=0d0 
           do i=1,ii
           k=j+(i-1)*nc
           dou(itau,k)=0d0 
           x=xx(i) 
           spec=0d0 
       if(dmu.le.dmu1.and.dmu.gt.dmu2) then 
          up=rr 
          down=1d0-t 
          diff=dabs(up-down)
          if(diff.gt.eps)
     &    call RADINT(SPEC,X,TAU,TAUT,T,RR,DMU,UP,DOWN,TIN,0) 
           dou(itau,k)=spec
       elseif(dmu.le.dmu2.and.dmu.gt.dmu3) then
           up=t-1d0
           down=0d0 
          diff=dabs(up-down)
          if(diff.gt.eps)
     &    call RADINT(SPEC,X,TAU,TAUT,T,RR,DMU,UP,DOWN,TIN,1) 
           dou(itau,k)=spec
           up=rr
           down=t-1d0 
          diff=dabs(up-down)
          if(diff.gt.eps)
     &    call RADINT(SPEC,X,TAU,TAUT,T,RR,DMU,UP,DOWN,TIN,0) 
           dou(itau,k)=dou(itau,k)+spec
       elseif(dmu.le.dmu3) then
           up=rr 
           down=0d0 
          diff=dabs(up-down)
c          if(diff.gt.eps)
          call RADINT(SPEC,X,TAU,TAUT,T,RR,DMU,UP,DOWN,TIN,1)
           dou(itau,k)=spec
       endif 
       sum2=sum2+dou(itau,k)
c end of x loop         
           enddo
       sum1=sum1+sum2*aang(j)*uang(j) 

c end of mu loop
         enddo 
       sum=sum+sum1*rr

c end of tau loop
       enddo
      return 
      END 

      SUBROUTINE RADINT(SPEC,X,TAU,TAUT,T,RR,DMU,UP,DOWN,TIN,IKEY) 
      IMPLICIT REAL*8(A-H,O-Z)
      PARAMETER(NGO=3) 
      COMMON/QQWGO/UGO(NGO),AGO(NGO) 
      COMMON/QWITCH/IGROUND,ICHAND,KEYRAY,IREDF,IDISTF,IGEOM,IPOL,ISP
      DATA PI/3.141592653589792384D0/,CMIN/7D1/,EPS/1D-10/   

      t2=t*t 
      dms=dsqrt(1d0-dmu**2) 

      if(rr.lt.eps) then 
      ngr=1 
      else 
      ngr=ngo 
      endif 
c integral over r  
      sum1=0d0 
         do kr=1,ngr 
         r=down+ugo(kr)*(up-down) 
         r2=r*r 
         rt=r*t*2d0 
         if(IKEY.eq.0) then
         phimax=dacos((1-r2-t2)/rt) 
         else 
         phimax=pi
         endif 

         sumphi=0d0 
c integral over phi
           do kph=1,ngo 
           phi=phimax*ugo(kph) 
           dis2=r2+t2+rt*dcos(phi) 
           dis=dsqrt(dis2) 
           cosph0=(dis2+t2-r2)/(2d0*dis*t) 
           if(IGEOM.eq.1) taupr=dmu/dms*
     &      (dis*cosph0-dsqrt(1d0-dis2+dis2*cosph0**2))
           if(IGEOM.eq.2) taupr=dmu*
     &      (dis*dms*cosph0-dsqrt(1d0-dis2+dis2*dms**2*cosph0**2))
           call TEMPH(TIN,DIS,TDIS)   
           blint=PLANCK(X,TDIS) 
           ab=0d0
           detmu=tau*(taut-taupr)/dmu
           if(detmu.lt.cmin) ab=dexp(-detmu) 
           sumphi=sumphi+ago(kph)*blint*ab 
           enddo 
           sumphi=sumphi*phimax 
           sum1=sum1+ago(kr)*r*sumphi 
         enddo 
         if(ngr.ne.1) then 
         sum1=sum1*(up-down)*2d0/rr**2
         else 
         sum1=sumphi
         endif 
         
      spec=sum1/pi  
      RETURN
      END 

      REAL FUNCTION PLANCK*8(X,TDIS) 
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/QWPH/YPH,YDISK 
      DATA EPS/1D-10/ 
      YD=511D0/TDIS 
      YX=X*YD 
      EYPH=0D0
      PLANCK=0D0 
      IF(YX.LT.7D1.and.YX.gt.1D-1) then 
      EYPH=DEXP(-YX)
      PLANCK=X**3*EYPH/(1D0-EYPH)
      elseif(YX.le.1D-1) then 
      SUM=0D0 
      ACC=1D0 
      DN=1D0 
      DO WHILE(ACC.gt.EPS)  
      SUM=SUM+ACC
      DN=DN+1D0 
      ACC=ACC*YX/DN 
      ENDDO
      PLANCK=X**2/(SUM*YD) 
      endif 
      PLANCK=PLANCK*YPH**3 
      RETURN
      END

      SUBROUTINE TEMPH(TIN,DIS,TDIS)
      IMPLICIT REAL*8(A-H,O-Z)
      DATA ALPHA/0.75D0/    
      TDIS=TIN*DIS**(-ALPHA)
      RETURN
      END

c multi-colour disk 
      SUBROUTINE BBMULTI(ENER,DIMULT,NN) 
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/WAGUE/UL(40),AL(40),AU(40),AE(40),N0
      COMMON/QWPH/YPH,YDISK 
      DIMENSION ENER(1),DIMULT(1)
      DATA EPS/1D-10/,CMAX/7D1/
      DATA C83/2.6666666667D0/,C53/1.666666667D0/,C13/0.333333333D0/    
      CONST=C83
      DO 10 IX=1,NN 
      X=ENER(IX) 
      YX=X*YPH 
      EYT=0D0
      IF(YX.lt.CMAX) EYT=DEXP(-YX) 
      SUM=0D0 
      DO 20 IA=1,N0  
      TT=UL(IA)+YX  
      EY=1D0 
      IF(TT.lt.CMAX) EY=1D0/(1D0-DEXP(-TT))    
      SUM=SUM+AL(IA)*EY*TT**C53 
   20 CONTINUE 
      DIMULT(IX)=SUM*CONST*YX**C13*EYT  
   10 CONTINUE 
      RETURN
      END

c***********************************************************************
C Incident black body spectrum:
      SUBROUTINE BLACKB(ENER,BLINT,N) 
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/QWPH/YPH,YDISK 
      DIMENSION ENER(1),BLINT(1) 
      DO 10 IX=1,N 
      X=ENER(IX)
      EYPH=0D0
      YX=YPH*X 
      IF(YX.LT.7D1) EYPH=DEXP(-YX) 
      BLINT(IX)=YX**3*EYPH/(1D0-EYPH)
   10 CONTINUE 
      RETURN
      END 
 
c*********************************************************************
c     calculates Compton reflected spectrum 
c     using Green's functions of Magdziarz & Zdziarski 1995 
c     and angle averaged incoming radiation 
c*********************************************************************
      subroutine MZrefl(x,compsp,xr,refl,angle,
     + par_refl,nrefl,jmax,jref,nc,kk,ioniz,irefl) 
c     compsp - angle averaged comptonized spectrum (flux)
c     xr - energies (1:jref) where refl should be computed
c     refl  - reflected spectrum (intensities) at different angles
c     angle - vector of angles (nc)
      PARAMETER(nem=5000) 
      implicit real*8 (a-h,o-z)
      integer jmax,jref,nc,ii,ia,nmin,nmax,i,kk,k,ioniz,irefl
      real*4 spref(nem),xref(nem),spco
      real*4 gamma,ab_,ab0,temp,xi,xm,ec
      real*8 x(jmax),xr(jref),compsp(jmax),refl(kk),angle(nc)
      real*4 cmin,par_refl(nrefl)
      common/wcompsp/XLOG(nem),DSPXX(nem),F2(nem),ii
      data cmin/1e-30/
c alog(1 keV/511) alog(10/511)
      data XMIN/-6.236/,XMAX/-3.934/,ec/0./
      external spco

c Fe abudance
      ab0=par_refl(1)
c  ab_       - abundance of elements heavier than He relative to
c              the solar abundances
      ab_=par_refl(2)
c  xi        - disk ionization parameter = L/nR^2
      xi=par_refl(3)
c  temp      - disk temperature (for ionization balance, Tmax=10^6 K)
      temp=par_refl(4)

      do i=1,jref
      xref(i)=sngl(xr(i))
      enddo
 
c put here calculation of spline for input spectrum
      ii=jmax
      do 212 i=1,jmax
      dspxx(i)=dlog(compsp(i)+cmin)
      xlog(i)=dlog(x(i))
  212 continue
      CALL QSPLINE(XLOG,DSPXX,JMAX,1d30,1d30,F2)

      if(ioniz.eq.0) then 
      gamma=2. 
      else 
c compute here photon spectral slope gamma in interval 1-10 keV
      CALL DLOCATE(XLOG,JMAX,XMIN,NMIN)
      NMIN=NMIN+1
      CALL DLOCATE(XLOG,JMAX,XMAX,NMAX)
      IF(NMIN.ge.NMAX) NMIN=NMAX-1
      CALL DLININT(DSPXX,XLOG,A,B1,NMIN,NMAX,JMAX)
      gamma=1.-sngl(A)
      endif
 
      do ia=1,nc
      xm=sngl(angle(ia))
c call compton reflection for a given angle
c spref gives  x F_x
c static reflection 
      if(irefl.eq.0) then 
      if(ioniz.eq.1) 
     &  call greeion(xref,spref,xm,ab_,ab0,spco,gamma,ec,temp,xi,jref)
      if(ioniz.eq.0) 
     &  call greeneu(xref,spref,xm,ab_,ab0,spco,gamma,ec,jref)
c relativistic disk reflection 
      else 
      call relrefl(xref,spref,xm,gamma,spco, 
     + par_refl,nrefl,jref,ioniz) 
      endif 

       do ix=1,jref
c intensity = flux/cos_angle
       k=ia+(ix-1)*nc
       refl(k)=dble(spref(ix))/xr(ix)/angle(ia)
       enddo
      enddo
 
      return
      end


c***********************************************************************
c compute n(1/lambda)*lambda^(-2) = F(1/lambda)/lambda
      real function spco(y,gamma,ec) 
      implicit real*8 (a-h,o-z)
      integer nem 
      real y,gamma,ec 
      PARAMETER(nem=5000) 
c      real XLOG(nem),DSPXX(nem),F2(nem)
      common/wcompsp/XLOG(nem),DSPXX(nem),F2(nem),II 
c use of the spline computed 
      CALL QSPLINT(XLOG,DSPXX,F2,II,dlog(dble(1./y)),YYY)
      spco=exp(sngl(yyy))/y 
      return
      end

c**************************************************
C linear interpolation:
C abciss:  XX
C ordinate: FUNLOG
C A,B: coefficients of the linear function Ax+B
      SUBROUTINE DLININT(FUNLOG,XLOG,A,B,NMIN,NMAX,II)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION FUNLOG(II),XLOG(II)
      DATA ERR/1D30/
      IF(NMAX.le.NMIN) THEN
      A=ERR
      B=ERR
      ELSE
      S=NMAX-NMIN+1
      SX=0. 
      SY=0. 
      SXX=0. 
      SXY=0. 
           DO 15 IX=NMIN,NMAX
           X=XLOG(IX)
           DI=FUNLOG(IX)
           SX=SX+X
           SY=SY+DI
           SXX=SXX+X*X
           SXY=SXY+DI*X
   15      CONTINUE
      DELTA=S*SXX-SX*SX
      B=(SXX*SY-SX*SXY)/DELTA
      A=(S*SXY-SX*SY)/DELTA
      ENDIF
      RETURN
      END
c***********************************************************************
      SUBROUTINE QSPLINE(X,Y,N,YP1,YPN,Y2)
      IMPLICIT REAL*8(A-H,O-Z) 
      PARAMETER (NMAX=500)
      DIMENSION X(N),Y(N),Y2(N),U(NMAX)
      IF (YP1.GT..99D30) THEN
        Y2(1)=0.
        U(1)=0.
      ELSE
        Y2(1)=-0.5
        U(1)=(3./(X(2)-X(1)))*((Y(2)-Y(1))/(X(2)-X(1))-YP1)
      ENDIF
      DO 11 I=2,N-1
        SIG=(X(I)-X(I-1))/(X(I+1)-X(I-1))
        P=SIG*Y2(I-1)+2.
        Y2(I)=(SIG-1.)/P
        U(I)=(6.*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1))
     *      /(X(I)-X(I-1)))/(X(I+1)-X(I-1))-SIG*U(I-1))/P
11    CONTINUE
      IF (YPN.GT..99D30) THEN
        QN=0.
        UN=0.
      ELSE
        QN=0.5
        UN=(3./(X(N)-X(N-1)))*(YPN-(Y(N)-Y(N-1))/(X(N)-X(N-1)))
      ENDIF
      Y2(N)=(UN-QN*U(N-1))/(QN*Y2(N-1)+1.)
      DO 12 K=N-1,1,-1
        Y2(K)=Y2(K)*Y2(K+1)+U(K)
12    CONTINUE
      RETURN
      END


      SUBROUTINE QSPLINT(XA,YA,Y2A,N,X,Y)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION XA(N),YA(N),Y2A(N)
      KLO=1
      KHI=N
1     IF (KHI-KLO.GT.1) THEN
        K=(KHI+KLO)/2
        IF(XA(K).GT.X)THEN
          KHI=K
        ELSE
          KLO=K
        ENDIF
      GOTO 1
      ENDIF
      H=XA(KHI)-XA(KLO)
      IF (H.EQ.0.) PAUSE 'Bad XA input.'
      A=(XA(KHI)-X)/H
      B=(X-XA(KLO))/H
      Y=A*YA(KLO)+B*YA(KHI)+
     *      ((A**3-A)*Y2A(KLO)+(B**3-B)*Y2A(KHI))*(H**2)/6.
      RETURN
      END
c*********************************************
c polinomial interpolation and extrapolation    
      SUBROUTINE POLINTQ(XA,YA,N,X,Y,DY)
      IMPLICIT REAL*8(A-H,O-Z) 
      PARAMETER (NMAX=10) 
      DIMENSION XA(N),YA(N),C(NMAX),D(NMAX)
      NS=1
      DIF=ABS(X-XA(1))
      DO 11 I=1,N 
        DIFT=ABS(X-XA(I))
        IF (DIFT.LT.DIF) THEN
          NS=I
          DIF=DIFT
        ENDIF
        C(I)=YA(I)
        D(I)=YA(I)
11    CONTINUE
      Y=YA(NS)
      NS=NS-1
      DO 13 M=1,N-1
        DO 12 I=1,N-M
          HO=XA(I)-X
          HP=XA(I+M)-X
          W=C(I+1)-D(I)
          DEN=HO-HP
          IF(DEN.EQ.0.D0)PAUSE 'In POLINTQ'
          DEN=W/DEN
          D(I)=HP*DEN
          C(I)=HO*DEN
12      CONTINUE
        IF (2*NS.LT.N-M)THEN
          DY=C(NS+1)
        ELSE
          DY=D(NS)
          NS=NS-1
        ENDIF
        Y=Y+DY
13    CONTINUE
      RETURN
      END

c**************************************************
c linear interpolation 
      SUBROUTINE DLINEAR(X1,Y1,X2,Y2,X,Y) 
      IMPLICIT REAL*8(A-H,O-Z) 
      DX=X2-X1 
      IF(DX.eq.0D0) PAUSE 'In DLINEAR'
      T=(Y2-Y1)/DX 
      Y=Y1+T*(X-X1) 
      RETURN
      END

c**********************************************************************
      SUBROUTINE QHUNT(XX,N,X,JLO)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION XX(N)
      LOGICAL ASCND
      ASCND=XX(N).GT.XX(1)
      IF(JLO.LE.0.OR.JLO.GT.N)THEN
        JLO=0
        JHI=N+1
        GO TO 3
      ENDIF
      INC=1
      IF(X.GE.XX(JLO).EQV.ASCND)THEN
1       JHI=JLO+INC
        IF(JHI.GT.N)THEN
          JHI=N+1
        ELSE IF(X.GE.XX(JHI).EQV.ASCND)THEN
          JLO=JHI
          INC=INC+INC
          GO TO 1
        ENDIF
      ELSE
        JHI=JLO
2       JLO=JHI-INC
        IF(JLO.LT.1)THEN
          JLO=0
        ELSE IF(X.LT.XX(JLO).EQV.ASCND)THEN
          JHI=JLO
          INC=INC+INC
          GO TO 2
        ENDIF
      ENDIF
3     IF(JHI-JLO.EQ.1)RETURN
      JM=(JHI+JLO)/2
      IF(X.GT.XX(JM).EQV.ASCND)THEN
        JLO=JM
      ELSE
        JHI=JM
      ENDIF
      GO TO 3
      END
c**************************************************
      SUBROUTINE DLOCATE(XX,N,X,J)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION XX(N)
      JL=0
      JU=N+1
10    IF(JU-JL.GT.1)THEN
        JM=(JU+JL)/2
        IF((XX(N).GT.XX(1)).EQV.(X.GT.XX(JM)))THEN
          JL=JM
        ELSE
          JU=JM
        ENDIF
      GO TO 10
      ENDIF
      J=JL
      RETURN
      END
c**************************************************

      SUBROUTINE AGUEGALA(NL,E)             
      IMPLICIT REAL*8(A-H,O-Z)          
      COMMON /WAGUE/ UL(40),AL(40),AU(40),AE(40),N0
      DLE=-DLOG(E)                           
      CALL QFQLAG0(NL,E,UL,AL)               
      DO 200 L=1,NL                          
C         IF(DABS(AL(L)).LT.E) RETURN        
         A=DEXP(-UL(L))                      
         AU(L)=A                             
         AE(L)=DEXP(-A*A)*2D0                
C         N0=L                                
  200 CONTINUE                               
       N0=NL
      RETURN                                 
      END                                    


      SUBROUTINE QDRGSDO(Z,W,NUM)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION Z(1),W(1)
      DATA ITRMAX/10/
      N=NUM
      IF (N.GT.1) GO TO 20
      IF (N.EQ.1) GO TO 10
      N=5
      GO TO 20
   10 Z(1)=0.D0
      W(1)=2.D0
      RETURN
   20 IND=MOD(N,2)
      K=N/2
      IF(IND.EQ.0) GOTO 40
      HP=1.D0
      FJ=1.D0
      DO 30 J=3,N,2
      FJ=FJ+2.D0
   30 HP=(FJ-1.D0)*HP/FJ
      HP2=HP*HP
      KIND=K+IND
      Z(KIND)=0.D0
      W(KIND)=HP2+HP2
   40 FN=DFLOAT(N)
      NP=N+1
      DO 100 I=1,K
      M=NP-I
      IF(I.EQ.1) X=1.D0-3.D0/(FN+2.D0)/FN
      IF(I.EQ.2) X=(X-1.D0)*5.25D0+1.D0
      IF(I.EQ.3) X=(X-Z(N))*1.6D0+X
      IF(I.GT .3) X=(X-Z(M+2))*3.D0+Z(M+3)
      ITER=0
   50 Y=X-1.D0
      ITER=ITER+1
      TY=2.D0*Y
      FJ=1.D0
      DP=Y
      P=X
      Q=1.D0
      DQ=1.D-0
      DO 60 J=2,N
      DQ=TY*Q+DQ+P
      Q=DQ+Q
      FJ=FJ+1.D0
      DP=(2.D0*FJ-1.D0)*Y*P+DP
   60 P=DP/FJ+P
      DIV=P/Q
      X=X-DIV
      FR=DIV/X
      IF(ITER.GT.ITRMAX) GOTO 190
      IF(DABS(FR).GT.1.D-16) GOTO 50
  190 CONTINUE
      Z(M)=X
      W(M)=2.D0/(1.D0-X*X)/Q/Q
      Z(I)=-X
      W(I)=W(M)
  100 CONTINUE
      RETURN
      END


C     SUBROUTINE FQLAG0(K,EPS,T,A) CALCULATE WEIGHTS AND POINTS 
C     K-POINTS GAUSSIAN QUADRATURE  
C     FROM INTERVAL  [0,INFINITY) WITH WEIGHT  EXP(-X).
C     EPS- RELATIVE ERROR OF POINTS CALCULATION 
C     T(1:K),A(1:K)-VECTORS OF POINTS AND WEIGHTS 
C     LIT.: A.K.PONOMARENKO,    LENINGRAD UNIVERSITY,1974 
C 
      SUBROUTINE QFQLAG0 (K,EPS,T,A)    
      IMPLICIT REAL*8(A-H,O-Z)         
      DIMENSION T(K),A(K)              
      I=0                              
      DP=K                             
      X=3.0D0/(1.0D0 + 2.4D0 * DP)     
      GO TO 30                         
   10 X=T(I) + 6.0D0 / (0.4D0 + DP)    
      GO TO 30                         
   20 J=I-1                            
      DJ=J                             
      X=T(I) + (T(I) - T(J)) * (1.0D0 + 2.55D0 * DJ)/1.9D0/DJ
   30 E2=1.0D0                         
      E3=1.0D0 - X                         
      IF (K.EQ.1) GO TO 41                 
      DO 40 J=2,K                          
       E1 = E2                             
      E2 = E3                              
      DJ = J                               
      DJ1= DJ - 1.0D0                      
   40 E3 = ((DJ + DJ1 - X)*E2 - DJ1 * E1)/DJ
   41 E1=DP*(E3-E2)                         
      E2 = E3/E1                            
      X  = X*(1.0D0 - E2)                   
      IF(DABS(E2)-EPS) 50,30,30             
   50 I=I+1                                 
      T(I)=X                                
      A(I)= X/E1/E1                         
      IF (I-1) 55,10,60                      
   55 continue 
      RETURN                                
   60 IF (I-K) 20,70,55                     
   70 RETURN                                
      END                                   


C SUBROUTINE  DBESK (X,N,BK,X1,NOUT) CALCULATE 
C FUNCTION CONNECTED WITH K(N,X) -- BESSEL FUNCTIONS OF ORDER 
C FROM  0 TO N, WHICH PUT TO THE FIRST N+1 ELEMENTS OF MASSIVE BK.
C C   X1 - BOUNDARY OF DIFFERENT  TYPES OF CALCULATIONS
C        (5<=X1<=8). 
C IF  0<=X<=X1  IS CALCULATED  K(N,X), IF   X>X1
C IS CALCULATED FUNCTION:
C
C    K(N,X)*EXP(X)*SQRT(2*X/PI), WHERE PI=3.14... 
C
C USE SUBROUTINES:
C
C    DBI0LE(X)      -8<=X<=8  
C    DBI1LO(X)      -8<=X<=8      DBK0GS(X)    X>=5 
C    DBK0LE(X)       0< X<=8      DBK1GS(X)    X>=5 
C    DBK1LO(X)       0< X<=8  
C
      SUBROUTINE QDBESK(X,N,BK,X1,NOUT)        
      IMPLICIT REAL*8(A-H,O-Z)                
      DIMENSION BK(N) 
      DATA EILER /0.57721 56649 01532 86061D0/
      DATA C1 /1D0/                           
C CONTROL OF ORDER:
      IF (N.GE.1) GO TO 10                    
C REACTION ON THE NEGATIVE ORDER:
                              stop 
c                              RETURN           
C CONTROL OF NEGATIVE ARGUMENT (X):
   10 IF (X) 13,15,20              
C REACTION ON THE NEGATIVE ARGUMENT 
   13                         stop 
c                              RETURN           
C CALCULATION IF X=0
   15 BK(1)=1D0     
      DO 17 I=2,N   
   17   BK(I)=0D0 
      RETURN        
C       CALCULATION OF K(0,X)
   20 IF (X.LE.X1) BK(1)=DBK0LE(X) 
      IF (X.LE.X1) BK(1)=BK(1)-(EILER+DLOG(X/2D0))*DBI0LE(X)
      IF (X.GT.X1) BK(1)=DBK0GS(X)                          
      IF (N.EQ.1) RETURN                                    
C CALCULATION OF K(1,X)
      IF (X.LE.X1) BK(2)=DBK1LO(X) 
      IF (X.LE.X1) ED=EILER+DLOG(X/2D0) 
      IF (X.LE.X1) BK(2)=ED*DBI1LO(X)+C1/X-BK(2) 
      IF (X.GT.X1) BK(2)=DBK1GS(X)               
      IF (N.EQ.2) RETURN                         
C CALCULATION OF K(N,X)
      DO 30 I=3,N      
        BK(I)=BK(I-2)+2*(I-2)*BK(I-1)/X  
   30 CONTINUE                             
      RETURN                               
      END                                        


C SUBROUTINE   DCHEB0 (A,X,N,KEY,NOUT)  CALCULATE 
C THE VALUES OF COEFFICIENTS OF EXPANSION OVER THE 
C CHEBYSHEV POLYNOMS
C IF ARE KNOWN FIRST  N+1 ELEMENTS OF VESTOR  A 
C 
C           F(N,X)=[K=0,N] A(K) T(L(KEY),X),     
C                                                
C IF        L        KEY                         
C         ______________                         
C                                                
C           N         0                          
C         2N+1        1                                    
C          2N         2                           
C           N (DISPL) 3,                          
C                                                 
C ALGORITHM KLENSHOW
C NOUT - NUMBER OF INFORMATION CHANAL OF  DCHEB0. 
C LITERATURE:Y.LUKE "SPECIAL MATHEMATICAL FUNCTIONS 
C                   AND ITS APPROXIMATIONS"
C                "MIR",MOSCOW  1980 608 P.
C  P.511
      FUNCTION QDCHEB0(A,X,N,KEY,NOUT)   
      IMPLICIT REAL*8 (A-H,O-Z)         
      DIMENSION A(*) 
      DATA CMAX /1D76/                  
      DATA C0,C2,C4 /0D0,2D0,4D0/       
C
C
C CONTROL OF NUMBER OF ITEMS
C
      IER=0  
      IF (N) 3,5,7
C        REACTION ON THE NEGATIVE NUMBER
    3                   QDCHEB0=CMAX    
                         IER=1           
    5 IF ((KEY.GE.0).AND.(KEY.LE.3)) GO TO 7       
C   CONTROL OF KEY
C  REACTION ON THE NONPOSITIVE VALUE OF KEY
                           QDCHEB0=1D76     
                           IER=1           
    7 IF (IER.NE.0) RETURN                 
      N1=N+1                               
      IF (KEY.EQ.0) Z=C2*X                 
      IF (KEY.EQ.1) Z=C4*X*X-C2            
      IF (KEY.EQ.2) Z=C4*X*X-C2            
      IF (KEY.EQ.3) Z=C4*X-C2              
      B0=C0                                
      B1=C0                                
      DO 10 I=1,N1                         
        K=N1-I                             
        B2=B1                              
        B1=B0                              
        B0=Z*B1-B2+A(K+1)                  
   10 CONTINUE                             
      IF (KEY.EQ.0) QDCHEB0=B0-X*B1         
      IF (KEY.EQ.1) QDCHEB0=X*(B0-B1)       
      IF (KEY.EQ.2) QDCHEB0=B0-B1*Z/C2      
      IF (KEY.EQ.3) QDCHEB0=B0-B1*Z/C2      
      RETURN                               
      END                                       



C FUNCTION  DBK0LE(X) CALCULATE 
C BESSEL FUNCTION   K(0,X) + I(0,X)*[ EILER+DLOG(X/2) ] 
C IF  X IS FROM INTERVAL [0,+8]
C  ( EILER=0.57721 56649 01532 86061 ).
C THROUGH EXPANSIONS OVER THE CHEBYSHEV POLYNOMS OF EVEN ORDER
C LITERATURE:Y.LUKE "SPECIAL MATHEMATICAL FUNCTIONS 
C                   AND ITS APPROXIMATIONS"
C                "MIR",MOSCOW  1980 608 P.
C TABLE  9.5  ; P.350
C COEFFICIENTS C
      FUNCTION DBK0LE(X)    
      IMPLICIT REAL*8 (A-H,O-Z)   
      DIMENSION B(19)             
      DATA NOUT /6/               
      DATA B / 240.27705 96407 20389 10102D0, 
     1         369.47407 39728 67282 63764D0, 
     2         169.97341 16984 01148 04378D0, 
     3          49.02046 37772 63439 39371D0, 
     4           9.38849 73252 68442 32756D0, 
     5           1.25947 97636 67703 58618D0, 
     6           0.12377 69641 14924 54118D0, 
     7           0.00924 43098 62866 90621D0, 
     8           0.00054 06238 96492 55807D0, 
     9           0.00002 53737 96028 08704D0, 
     *           0.00000 09754 78302 83898D0, 
     1           0.00000 00312 49571 77932D0, 
     2           0.00000 00008 46434 70610D0, 
     3           0.00000 00000 19628 88451D0, 
     4           0.00000 00000 00393 96098D0, 
     5           0.00000 00000 00006 90835D0, 
     6           0.10673D-15, 0.146D-17, 0.2D-19  /
C                                                  
      W=X*0.125D0                                  
      DBK0LE=QDCHEB0(B,W,18,2,NOUT)                 
      RETURN                                       
C                                                  
      END                                          



C FUNCTION  DBK1LO(X) CALCULATE 
C BESSEL FUNCTION 
C            1/X + I(1,X)*[EILER+LN(X/2)]- K(1,X)       
C  IF X IS FROM INTERVAL [0,+8]
C ( EILER=0.57721 56649 01532 86061 ).
C THROUGH EXPANSIONS OVER THE CHEBYSHEV POLYNOMS OF ODD  ORDER
C LITERATURE:Y.LUKE "SPECIAL MATHEMATICAL FUNCTIONS 
C                   AND ITS APPROXIMATIONS"
C                "MIR",MOSCOW  1980 608 P.
C TABLE  9.6  ; P.352
C COEFFICIENTS C 
      FUNCTION DBK1LO(X)                      
      IMPLICIT REAL*8 (A-H,O-Z)               
      DIMENSION C(18)                         
      DATA NOUT /6/                           
      DATA C / 418.88944 61663 96890 97522D0, 
     1         249.89554 90428 68080 38961D0, 
     2          91.18031 93387 41787 75763D0, 
     3          21.44499 50539 62240 43921D0, 
     4           3.43841 53928 80464 59793D0, 
     5           0.39484 60929 40938 23432D0, 
     6           0.03382 87455 26884 19281D0, 
     7           0.00223 57203 34170 88760D0, 
     8           0.00011 71310 22460 84561D0, 
     9           0.00000 49754 27122 13645D0, 
     *           0.00000 01746 04931 76984D0, 
     1           0.00000 00051 43294 11806D0, 
     2           0.00000 00001 28903 39664D0, 
     3           0.00000 00000 02780 94119D0, 
     4           0.00000 00000 00052 17097D0, 
     5           0.00000 00000 00000 85869D0, 
     6           0.00000 00000 00000 01250D0, 
     7           0.00000 00000 00000 00016D0 /
      W=X*0.125D0                             
      DBK1LO=QDCHEB0(C,W,17,1,NOUT)            
      RETURN                                  
C                                             
      END                                     



C FUNCTION  DBK0GS(X) CALCULATE 
C BESSEL FUNCTION K(0,X)*SQRT(2*X/PI)*EXP(X) IF X>=5
C THROUGH EXPANSIONS OVER THE SHIFTED(DISPLACED) CHEBYSHEV POLYNOMS 
C LITERATURE:Y.LUKE "SPECIAL MATHEMATICAL FUNCTIONS 
C                   AND ITS APPROXIMATIONS"
C                "MIR",MOSCOW  1980 608 P.
C TABLE  9.5  ; P.351
C COEFFICIENTS D
      FUNCTION DBK0GS(X)      
      IMPLICIT REAL*8 (A-H,O-Z)  
      DIMENSION D(21)            
      DATA NOUT /6/              
      DATA D / 0.98840 81742 30825 80035D0, 
     1        -0.01131 05046 46928 28069D0, 
     2         0.00026 95326 12762 72369D0, 
     3        -0.00001 11066 85196 66535D0, 
     4         0.00000 06325 75108 50049D0, 
     5        -0.00000 00450 47337 64110D0, 
     6         0.00000 00037 92996 45568D0, 
     7        -0.00000 00003 64547 17921D0, 
     8         0.00000 00000 39043 75576D0, 
     9        -0.00000 00000 04579 93622D0, 
     *         0.00000 00000 00580 81063D0, 
     1        -0.00000 00000 00078 83236D0, 
     2         0.00000 00000 00011 36042D0, 
     3        -0.00000 00000 00001 72697D0, 
     4  0.27545D-15,   -0.4589D-16,    0.796D-17,  -0.143D-17,
     5  0.27D-18,      -0.5D-19,       0.1D-19  /             
C
C
      W=5D0/X                                                 
      DBK0GS=QDCHEB0(D,W,20,3,NOUT)                           
      RETURN                                                 
C
      END                                                    



C FUNCTION  DBK1GS(X) CALCULATE 
C BESSEL FUNCTION K(1,X)*SQRT(2*X/PI)*EXP(X) IF X>=5
C THROUGH EXPANSIONS OVER THE SHIFTED(DISPLACED) CHEBYSHEV POLYNOMS 
C LITERATURE:Y.LUKE "SPECIAL MATHEMATICAL FUNCTIONS 
C                   AND ITS APPROXIMATIONS"
C                "MIR",MOSCOW  1980 608 P.
C TABLE  9.6  ; P.351
C     COEFFICIENTS   E  
      FUNCTION DBK1GS(X)                      
      IMPLICIT REAL*8 (A-H,O-Z)               
      DIMENSION D(21)                         
      DATA NOUT /6/                           
      DATA D / 1.03595 08587 72358 33071D0,   
     1         0.03546 52912 43331 11380D0,   
     2        -0.00046 84750 28166 88856D0,   
     3         0.00001 61850 63810 05343D0,   
     4        -0.00000 08451 72048 12368D0,   
     5         0.00000 00571 32218 10284D0,   
     6        -0.00000 00046 45554 60661D0,   
     7         0.00000 00004 35417 33857D0,   
     8        -0.00000 00000 45757 29704D0,   
     9         0.00000 00000 05288 13281D0,   
     *        -0.00000 00000 00662 61293D0,   
     1         0.00000 00000 00089 04792D0,   
     2        -0.00000 00000 00012 72607D0,   
     3         0.00000 00000 00001 92086D0,   
     4 -0.30450D-15,    0.5046D-16,   -0.871D-17,   0.156D-17, 
     5 -0.29D-18,       0.6D-19,      -0.1D-19  / 
C
C
      W=5D0/X                                     
      DBK1GS=QDCHEB0(D,W,20,3,NOUT)                
      RETURN                                      
C                                                 
      END                                         



C FUNCTION  DBI0LE(X) CALCULATE 
C BESSEL FUNCTION I(0,X) IF X IS FROM INTERVAL [-8,+8]
C THROUGH EXPANSIONS OVER THE CHEBYSHEV POLYNOMS OF EVEN ORDER
C LITERATURE:Y.LUKE "SPECIAL MATHEMATICAL FUNCTIONS 
C                   AND ITS APPROXIMATIONS"
C                "MIR",MOSCOW  1980 608 P.
C TABLE  9.5  ; P.350
      FUNCTION DBI0LE(X)                      
      IMPLICIT REAL*8 (A-H,O-Z)               
      DIMENSION A(19)                         
      DATA NOUT /6/                           
      DATA A / 127.73343 98121 81083 56301D0, 
     1         190.49432 01727 42844 19322D0, 
     2          82.48903 27440 24099 61321D0, 
     3          22.27481 92424 62230 87742D0, 
     4           4.01167 37601 79348 53351D0, 
     5           0.50949 33654 39982 87079D0, 
     6           0.04771 87487 98174 13524D0, 
     7           0.00341 63317 66012 34095D0, 
     8           0.00019 24693 59688 11366D0, 
     9           0.00000 87383 15496 62236D0, 
     *           0.00000 03260 91050 57896D0, 
     1           0.00000 00101 69726 72769D0, 
     2           0.00000 00002 68828 12895D0, 
     3           0.00000 00000 06096 89280D0, 
     4           0.00000 00000 00119 89083D0, 
     5           0.00000 00000 00002 06305D0, 
     6           0.3132D-16, 0.42D-18, 0.1D-19  /
C                                                
      W=X*0.125D0                                
      DBI0LE=QDCHEB0(A,W,18,2,NOUT)               
      RETURN                                     
C                                                
      END                                        



C FUNCTION  DBI1LO(X) CALCULATE 
C BESSEL FUNCTION I(1,X) IF X IS FROM INTERVAL [-8,+8]
C THROUGH EXPANSIONS OVER THE CHEBYSHEV POLYNOMS OF ODD  ORDER
C LITERATURE:Y.LUKE "SPECIAL MATHEMATICAL FUNCTIONS 
C                   AND ITS APPROXIMATIONS"
C                "MIR",MOSCOW  1980 608 P.
C TABLE  9.6  ; P.352
      FUNCTION DBI1LO(X)                      
      IMPLICIT REAL*8 (A-H,O-Z)               
      DIMENSION A(18)                         
      DATA NOUT /6/                           
      DATA A / 220.60142 69235 23778 56112D0, 
     1         125.35426 68371 52356 46451D0, 
     2          42.86523 40931 28256 85130D0, 
     3           9.45300 52294 34910 53517D0, 
     4           1.42965 77090 76213 46814D0, 
     5           0.15592 42954 76256 29116D0, 
     6           0.01276 80490 81733 88545D0, 
     7           0.00081 08879 00690 69214D0, 
     8           0.00004 10104 61938 23750D0, 
     9           0.00000 16880 42203 43687D0, 
     *           0.00000 00575 86950 54206D0, 
     1           0.00000 00016 53453 53976D0, 
     2           0.00000 00000 40484 76606D0, 
     3           0.00000 00000 00854 96289D0, 
     4           0.00000 00000 00015 72708D0, 
     5           0.00000 00000 00000 25419D0, 
     6           0.00000 00000 00000 00364D0, 
     7           0.00000 00000 00000 00005D0 /
      W=X*0.125D0                             
      DBI1LO=QDCHEB0(A,W,17,1,NOUT)  
      RETURN                                  
C                                             
      END                                     

C--------------------------------------------------------------------------
C--------------------------------------------------------------------------

