SUBROUTINE COMSET
C     =================
C
C     sets up necessary parameters for treating the Compton scattering
c
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'MODELQ.FOR'
      dimension freqi(mfreq)
      parameter (xcon=8.0935d-21,YCON=1.68638E-10)
      parameter (t15=1.d-15)
      common/auxcbc/cden1m(mdepth),cden10(mdepth),
     *              cden2m(mdepth),cden20(mdepth)
      common/comgfs/gfm(mfreq,mdeptc),gfp(mfreq,mdeptc)
      DIMENSION PL(MDEPTH),PLM(MDEPTH)
C
      if(icompt.le.0) go to 100
      nmuc=3
      nsti=0
      nedd=3
C
C     frequency-dependent universal parameters
C
      do ij=1,nfreq
         cder10(ij)=0.
         cder1p(ij)=0.
         cder1m(ij)=0.
         cder20(ij)=0.
         cder2p(ij)=0.
         cder2m(ij)=0.
         iji=nfreq-kij(ij)+1
         ijorig(iji)=ij
         freqi(iji)=freq(ij)
         fr=freqi(iji)
         bnus(iji)=two*xcon*fr/(bn*(fr*1.d-15)**3)
      end do
C
      ij=1
      dlnfr(ij)=log(freqi(ij+1)/freqi(ij))
      do ij=2,nfreq-1
         dlnfr(ij)=log(freqi(ij+1)/freqi(ij))
         delp=dlnfr(ij)
         delm=dlnfr(ij-1)
         del0=delp+delm
         cd0=two/del0
         cder2m(ij)=cd0/delm
         cder2p(ij)=cd0/delp
         cder20(ij)=-cder2m(ij)-cder2p(ij)
      end do
c
      do ij=1,nfreq-1
         frj0=freqi(ij)
         frjp=freqi(ij+1)
         frz=sqrt(frj0*frjp)
         do id=1,nd
C           to avoid over/underflow problems:
            IF(HK*FRJ0/TEMP(ID).LT.200.) THEN
              fjb0=un/(exp(hk*frj0/temp(id))-un)
             ELSE
              fjb0=0.
            ENDIF
            IF(HK*FRJP/TEMP(ID).LT.200.) THEN
               fjbp=un/(exp(hk*frjp/temp(id))-un)
             ELSE
              fjbp=0.
            ENDIF
          fjz0=fjb0*(bn*(frj0*t15)**3)
          fjzp=fjbp*(bn*(frjp*t15)**3)
          if(ichcoo.eq.0) then
             zj0=hk*frz/temp(id)
             dfjz=fjz0-fjzp
             dfjb=fjb0-fjbp
             fzz=un+fjbp-3./zj0
             aa=dfjz*dfjb
             bb=dfjz*fzz+fjzp*dfjb
             cc=fjzp*fzz-dfjz/dlnfr(ij)/zj0
           else
               e2=ycon*temp(id)
               zxxp=xcon*frjp*(un+fjbp)-3.*e2
               zxx0=xcon*frj0*(un+fjb0)-3.*e2
             dzxx=zxx0-zxxp
             dfjb=fjb0-fjbp
               dfjz=fjz0-fjzp
             aa=dfjz*dzxx
             bb=dfjz*zxxp+fjzp*dzxx
             cc=fjzp*zxxp-e2*dfjz/dlnfr(ij)
          end if
CXXX        to avoid division by zero:
          if(abs(aa).eq.0.and.abs(bb).eq.0.) then
              xx1=0.
            elseif(abs(aa).lt.1.e-7*abs(bb)) then
             xx1=-cc/bb
          else
             dd=bb*bb-4.*aa*cc
             if(dd.lt.0.) dd=0.
             dd=sqrt(dd)
             xx1=(dd-bb)*half/aa
             if(ichcoo.gt.0) then
                  xx2=-(dd+bb)*half/aa
                  dxx1=abs(xx1-half)
                  dxx2=abs(xx2-half)
                  if(dxx2.lt.dxx1) xx1=xx2
                  if((xx1.gt.1.).or.(xx1.lt.0.)) xx1=half
             end if
            end if
          delj(ij,id)=xx1
         end do
      end do
c
C     angle-dependent universal parameters
C
      call angset
c
C     frequency-dependent universal parameters
C
  100 continue
      do ij=1,nfreq
c
c  first-order expression
c
         if(knish.eq.0) then
            SIGEC(IJ)=SIGE*(un-two*freq(ij)*xcon)
c
C  Use full Klein-Nishina cross section (Rybicki & Lightman 1975):
c
          else
            xf=xcon*freq(ij)
            if(xf.lt.1.d-1) then
               SIGEC(IJ)=SIGE*(1.-xf*(2.-xf*(26./5.-xf*(13.3
     *                  -xf*(1144./35.-xf*(544./7.-xf*(3784./21.
     *                  -xf*(6148./15.-xf*(151552./165.
     *                  -xf*111872./55.)))))))))
              else if(xf.gt.1.d3) then
                 SIGEC(IJ)=SIGE*3./8./xf*(log(2.*xf)+0.5)
              else
                 SIGEC(IJ)=SIGE*0.75*((1.+xf)/xf**3*(2.*xf*(1.+xf)/
     *           (1.+2.*xf)-log(1.+2.*xf))+0.5*log(1.+2.*xf)/xf
     *           -(1.+3.*xf)/(1.+2.*xf)**2)
            endif
         endif
      end do
c
      if(icompt.le.0) return
      IJ=1
      IJO=ijorig(ij)
      DO ID=1,ND
         PLM(ID)=BNUE(IJO)/(EXP(HK/Temp(ID)*FREQ(IJO))-UN)
      END DO
C      
      DO IJ=2,NFREQ
       IJO=ijorig(ij)
         DO ID=1,ND
C           to avoid over/underflow problems:
            IF(HK/TEMP(ID)*FREQ(IJO).LT.200.) THEN
               PL(ID)=BNUE(IJO)/(EXP(HK/temp(ID)*FREQ(IJO))-UN)
            ELSE
               PL(ID)=PLM(ID)
            ENDIF
            PLM(ID)=PL(ID)
         END DO
      END DO
C      
      return
      end