C**********************************************************************
C*    SIMPLEX p lfc ym    ib          k  'SPECTS' 'PHIV'              *
C*                                ORIGINALLY CODED BY E.FENIMORE      *
C*  PBSG150.F.FORT77(MODEL) NI COMPTONIZED P-L O KUWAERU              *
C*        B.B. UNDER FOLLOW PROTECTION                                *
C*        OCULTED COMPTNIZE O IRERU                                   *
C**********************************************************************
      SUBROUTINE SPECTS(EINDX1,EINDX2,ITY,SCALE,E1,E2,ESTART,EEND,
     * EINDX3,EINDX4,PHOTON)
      IMPLICIT REAL*8 (A-H,O-Z)
C
C     THE PURPOSE OF SPECTS IS TO GIVE THE NUMBER OF PHOTON BETWEEN
C     E1 ENERGY AND E2.
C
      COMMON /PARMS/TS,TE,ODEPTH,EF,X
      COMMON /COEF/ COEFB(7,2),GAUSS(5,2),BREAK(7),CBACK(15,2)
      COMMON /NCOEF/ NUMGU,NUMBKG
      DATA NUMINT/1/,CFWSD/0.60056D0/
      DIMENSION COEFBU(7,2)
      PHOTON=0.
C
C     DO POWER LAW
C
      IF(ITY.EQ.1) THEN
         IF(ABS(EINDX1-1.).GE.(.0001)) THEN
            POWER=1.-EINDX1
            PLCNST=1./POWER
C-------------- KITAMOTO CHANGE --------
             IF(POWER.GT.250.)POWER=252.
            SUM1=E1**(POWER)
            SUM2=E2**(POWER)
            PHOTON=PLCNST*(SUM2-SUM1)
         ELSE
ce          PHOTON=ALOG(E2/E1)
            PHOTON=DLOG(E2/E1)
         ENDIF
C
C     DO LIANG THERMAL CYCLTRON
C
      ELSEIF(ITY.EQ.2) THEN
         DB=(E2-E1)/NUMINT
         DDB=DB/2.
         SUM=0.
         DO 211 I=1,NUMINT
         APBB=(2*I-1)*DDB+E1
         DO 211 II=1,NUMGU
         E=DDB*GAUSS(II,2)+APBB
         AE=4.5*E*EINDX1
         TERM=EXP(-AE**(1./3.))
         SUM=SUM+GAUSS(II,1)*TERM
211      CONTINUE
         PHOTON=SUM*DDB
C
C     DO THERMAL BREMSSTRAHLUNG (CLASSICAL)
C
      ELSEIF(ITY.EQ.3) THEN
         DB=(E2-E1)/NUMINT
         DDB=DB/2.
         SUM=0.
         DO 32 I=1,NUMINT
         APBB=(2*I-1)*DDB+E1
         DO 32 II=1,NUMGU
         E=DDB*GAUSS(II,2)+APBB
         AE=E/EINDX1
C---------  CHAPLIN STEVENS SPECTRUM  .... S.K. --------
           IF(EINDX2.GT.0.0001)THEN
            CALL CHASTE(E,EINDX1,EINDX2,CR)
            TERM=TERM*CR
           ENDIF
C--------------------------------------------------------
         TERM=EXP(-AE)/E
         SUM=SUM+GAUSS(II,1)*TERM
32       CONTINUE
         PHOTON=SUM*DDB
C
C     DO BLACK BODY LAW
C
      ELSEIF(ITY.EQ.4) THEN
C        BBCNST=EINDX1**(3.)
C        ARG1=E1/EINDX1
C        C11=EXP(-ARG1)
C        EXPQ1=1.
C        ARG=E2/EINDX1
C        C1=EXP(-ARG)
C        EXPQ=1.
C        SUM1=0.
C        SUM2=0.
C        DO 45 N=1,10
C        Q1=ARG1*N
C        Q21=Q1*Q1
C        POLY1=Q21+2.*Q1+2.
C        EXPQ1=C11*EXPQ1
C        SUM1=SUM1+EXPQ1*POLY1/(N**3)
C        Q=ARG*N
C        Q2=Q*Q
C        POLY=Q2+2.*Q+2.
C        EXPQ=C1*EXPQ
C        SUM2=SUM2+EXPQ*POLY/(N**3)
C45       CONTINUE
C         PHOTON=BBCNST*(SUM1-SUM2)
         DB=(E2-E1)/NUMINT
         DDB=DB/2.
         SUM=0.
         DO 46 I=1,NUMINT
         APBB=(2*I-1)*DDB+E1
         DO 46 II=1,NUMGU
         E=DDB*GAUSS(II,2)+APBB
         U=E/EINDX1
C----------- KITAMOTO CHANGE ------------
C        IF(U.GT.500.) U=500 --- X.LT.174.673 FOR EXP(X)
         IF(U.GT.170.) GO TO 46
         TERM=E*E/(EXP(U)-1.)
         SUM=SUM+GAUSS(II,1)*TERM
46       CONTINUE
         PHOTON=SUM*DDB
C
C     DO THERM BREM, MATTESON OR COMBINED MATTESON-JACOBSON
C
      ELSEIF(ITY.EQ.5.OR.ITY.EQ.6) THEN
         UNIT=EINDX1/100.
         PW1=-.333
         IF((1./EINDX1).GT..01)PW1=-.190
         PW2=-.31*UNIT**PW1
         PW3=.415
         IF((1./EINDX1).GT..01)PW3=.125
         COEFTB=1.04*UNIT**PW3
         IF(ITY.EQ.6) THEN
            DO 69 I=1,7
            COEFBU(I,1)=1.+COEFB(I,1)*EINDX1/511.
            COEFBU(I,2)=COEFB(I,2)/511.
69          CONTINUE
         ENDIF
         DB=(E2-E1)/NUMINT
         DDB=DB/2.
         SUM=0.
         DO 51 I=1,NUMINT
         APBB=(2*I-1)*DDB+E1
         DO 51 II=1,NUMGU
         E=DDB*GAUSS(II,2)+APBB
         AE=E/EINDX1
C------------ KITAMOTO ----------------
           IF(AE.GT.170.)AE=170.
         TERM=EXP(-AE)*AE**PW2/E
         IF(ITY.EQ.6)THEN
            IBST=7
            DO 60 WHILE(AE.GT.BREAK(IBST))
60          IBST=IBST-1
            TERM=TERM*(COEFBU(IBST,1)+COEFBU(IBST,2)*E)
         ENDIF
C---------  CHAPLIN STEVENS SPECTRUM  .... S.K. --------
           IF(EINDX2.GT.0.0001)THEN
            CALL CHASTE(E,EINDX1,EINDX2,CR)
            TERM=TERM*CR
           ENDIF
C--------------------------------------------------------
         SUM=SUM+GAUSS(II,1)*TERM
51       CONTINUE
         PHOTON=SUM*COEFTB*DDB
C
C     DO EXPONENTIAL
C
      ELSEIF(ITY.EQ.7) THEN
         SUM1=EXP(-E1/EINDX1)*EINDX1
         SUM2=EXP(-E2/EINDX1)*EINDX1
         PHOTON=SUM1-SUM2
C
C     DO SEGMENT OF A POWER LAW DEFINED BY TWO POINTS
C
      ELSEIF(ITY.EQ.8) THEN
ce       ARB1=ALOG(EINDX1)
ce       ARB2=ALOG(EINDX2)
         ARB1=DLOG(EINDX1)
         ARB2=DLOG(EINDX2)
ce       SLOPE=(ARB2-ARB1)/ALOG(EEND/ESTART)
ce       BB=ARB1-SLOPE*ALOG(ESTART)
         SLOPE=(ARB2-ARB1)/DLOG(EEND/ESTART)
         BB=ARB1-SLOPE*DLOG(ESTART)
         SLOPEP=SLOPE+1.
         IF(ABS(SLOPEP).LE.(.01)) THEN
ce          PHOTON=EXP(BB)*ALOG(E2/E1)
            PHOTON=EXP(BB)*DLOG(E2/E1)
         ELSE
            SUM=(E2**SLOPEP-E1**SLOPEP)/SLOPEP
            PHOTON=SUM*EXP(BB)
         ENDIF
C
C     DO GAUSSIAN LINE, NOTE THAT LOADFR ASSUME THAT GAUSSIAN IS
C     FORMULA NUM 9
C
      ELSEIF(ITY.EQ.9) THEN
         DB=(E2-E1)/NUMINT
         DDB=DB/2.
         SUM=0.
         W=CFWSD*EINDX2
         DO 102 I=1,NUMINT
         APBB=(2*I-1)*DDB+E1
         DO 102 II=1,NUMGU
         E=DDB*GAUSS(II,2)+APBB
         ZZ=(E-EINDX1)/W
         TERM=EXP(-ZZ**2)
         SUM=SUM+GAUSS(II,1)*TERM
102      CONTINUE
         PHOTON=SUM*DDB/(W*1.77245385)
C
C     SUMMING OF VARIOUS ABSORPTIONS IS NOT ALLOWED, IT ONLY IS USED
C     FOR MULTIPLICATION (IN PHIV)
C
      ELSEIF(ITY.GE.10.AND.ITY.LE.14) THEN
         PHOTON=0
C******* FOR FIREMAN ABSORP TION 84.1.31  S.K. *******
      ELSEIF(ITY.EQ.30) THEN
         PHOTON=0
C******* FOR USER ABSORPTION     84.1.31  S.K. *******
      ELSEIF(ITY.EQ.31) THEN
         PHOTON=0
C******* A*E**(-ALPH) ABSORPTION 84.1.31  S.K. *******
      ELSEIF(ITY.EQ.32) THEN
         PHOTON=0
C******* A*E**(-ALPH) ABSORPTION 84.1.31  S.K. *******
      ELSEIF(ITY.EQ.33) THEN
         PHOTON=0
C
C     DO EXP CUTOFF
C
      ELSEIF(ITY.EQ.15) THEN
         DB=(E2-E1)/NUMINT
         DDB=DB/2.
         SUM=0.
         DO 162 I=1,NUMINT
         APBB=(2*I-1)*DDB+E1
         DO 162 II=1,NUMGU
         E=DDB*GAUSS(II,2)+APBB
         TERM=EXP((EINDX1-E)/EINDX2)
         SUM=SUM+GAUSS(II,1)*TERM
162      CONTINUE
         PHOTON=SUM*DDB
C
C     DO DIFFUSE BACKGROUND, NOTE THAT FORSC SHOULD TO THE STER OF
C     THE SYSTEM
C
      ELSEIF(ITY.EQ.16) THEN
         DB=(E2-E1)/NUMINT
         DDB=DB/2.
         SUM=0.
         DO 172 I=1,NUMINT
         APBB=(2*I-1)*DDB+E1
         DO 172 II=1,NUMGU
         E=DDB*GAUSS(II,2)+APBB
ce       ENERA=ALOG10(E)
         ENERA=DLOG10(E)
         DO 171 K=2,NUMBKG
         IF(ENERA.LT.CBACK(K,2))GO TO 174
171      CONTINUE
         K=NUMBKG
174      K=K-1
         SLOPE=(CBACK(K,1)-CBACK(K+1,1))/(CBACK(K,2)-CBACK(K+1,2))
         TERM2=CBACK(K+1,1)+(ENERA-CBACK(K+1,2))*SLOPE
         TERM=10.**TERM2
         SUM=SUM+GAUSS(II,1)*TERM
172      CONTINUE
         PHOTON=SUM*DDB
      ELSEIF(ITY.EQ.35.OR.ITY.EQ.36) THEN
C******* COMPTONIZED B-B , PL*****  ******84.8.6    S.K./9.4 S.K.
C      EINDX1=KT(B.B) KEV(FOR 35)/ ALPH(P.L.) (FOR.36)
C      EINDX2=TAU     OPTICAL DEPTH
C      EINDX3=KT(ELECTRON) KEV
C
         PROB=1.-EXP(-EINDX2)
         GAMMA=1.0+7.8278E-3*EINDX3+4.5955E-5*EINDX3**2
         NPORI=20
         DB=E2-E1
         E=DB/2.+E1
         CALL BCOMP(EINDX1,GAMMA,PROB,E,NPORI,FOUT,ITY)
         PHOTON=FOUT*DB
      ELSEIF(ITY.EQ.38.OR.ITY.EQ.39) THEN
C******* COMPTONIZED B-B (OCULTATION)*****84.11.2    S.K./9.4 S.K.
C      EINDX1=KT(B.B) KEV(FOR 38)/ ALPH(P.L.) (FOR.39)
C      EINDX2=TAU     OPTICAL DEPTH
C      EINDX3=KT(ELECTRON) KEV
C
         PROB=1.-EXP(-EINDX2)
         GAMMA=1.0+7.8278E-3*EINDX3+4.5955E-5*EINDX3**2
         NPORI=20
         DB=E2-E1
         E=DB/2.+E1
         CALL BCOMPO(EINDX1,GAMMA,PROB,E,NPORI,FOUT,ITY)
         PHOTON=FOUT*DB
      ENDIF
C
C     DO RETURN
C
      PHOTON=PHOTON*SCALE
      RETURN
      END
      SUBROUTINE PHIV(EINDX1,EINDX2,ITY,SCALE,ENER,ESTART,EEND,
     &  EINDX3,EINDX4,SPHI)
      IMPLICIT REAL*8 (A-H,O-Z)
C
C     THE PURPOSE OF PHIV IS TO GIVE THE VALUE OF THE PHOTON SPECTRUM
C     AT ENER
C
      COMMON /PARMS/TS,TE,ODEPTH,EF,X
      COMMON /COEF/ COEFB(7,2),GAUSS(5,2),BREAK(7),CBACK(15,2)
      COMMON /NCOEF/ NUMGU,NUMBKG
      DATA NUMINT/1/,CFWSD/0.6D0/
      DIMENSION FELABS(4),FEKABS(4),COEFBU(7,2)
      DATA FELABS/13.6696,-2.39195,-.137648,0./
      DATA FEKABS/14.3456,-1.23491,-.418785,.0321662/
      DATA FEEDGE/7.111/
      DATA FE26ED,FE25ED/9.2,8.7/
ce    DATA CFWSD/0.6/
      DATA DIPOLE,PI,PISQ/2.8,3.141592653,1.77245385/
      IF(ITY.EQ.1) THEN
         SPHI=ENER**(-EINDX1)
      ELSEIF(ITY.EQ.2) THEN
         AE=4.5*ENER*EINDX1
         SPHI=EXP(-AE**.333333333)
      ELSEIF(ITY.EQ.3) THEN
         AE=ENER/EINDX1
         IF(AE.GT.500)AE=500
C----------- CHAPLINE STEVENS SPECTRUM..... S.K. --------
           CR=1.0000
           IF(EINDX2.GT.0.0001)CALL CHASTE(ENER,EINDX1,EINDX2,CR)
           SPHI=EXP(-AE)/ENER *CR
C--------------------------------------------------------
C        SPHI=EXP(-AE)/ENER
      ELSEIF(ITY.EQ.4) THEN
         U=ENER/EINDX1
         IF(U.GT.500)U=500
         SPHI=ENER*ENER/(EXP(U)-1.)
      ELSEIF(ITY.EQ.5.OR.ITY.EQ.6) THEN
         UNIT=.01*EINDX1
         PW1=-.3333
         IF((1./EINDX1).GT..01)PW1=-.190
         PW2=-.31*UNIT**PW1
         PW3=.415
         IF((1./EINDX1).GT..01)PW3=.125
         AE=ENER/EINDX1
         IF(AE.GT.500)AE=500
         COEFTB=1.04*UNIT**PW3
         SPHI=COEFTB*EXP(-AE)*AE**PW2/ENER
         IF(ITY.EQ.6) THEN
            DO 69 I=1,7
            COEFBU(I,1)=1.+COEFB(I,1)*EINDX1/511.
            COEFBU(I,2)=COEFB(I,2)/511.
69          CONTINUE
            IBST=7
            DO 61 WHILE(AE.GT.BREAK(IBST))
61          IBST=IBST-1
            SPHI=SPHI*(COEFBU(IBST,1)+COEFBU(IBST,2)*ENER)
         ENDIF
C----------- CHAPLINE STEVENS SPECTRUM..... S.K. --------
           CR=1.0000
           IF(EINDX2.GT.0.0001)CALL CHASTE(ENER,EINDX1,EINDX2,CR)
           SPHI=SPHI *CR
C--------------------------------------------------------
      ELSEIF(ITY.EQ.7) THEN
         AE=ENER/EINDX1
         IF(AE.GT.500)AE=500
         SPHI=EXP(-AE)
      ELSEIF(ITY.EQ.8) THEN
ce       ARB1=ALOG(EINDX1)
ce       ARB2=ALOG(EINDX2)
         ARB1=DLOG(EINDX1)
         ARB2=DLOG(EINDX2)
ce       SLOPE=(ARB2-ARB1)/ALOG(EEND/ESTART)
ce       B=ARB1-SLOPE*ALOG(ESTART)
ce       TERM=SLOPE*ALOG(ENER)+B
         SLOPE=(ARB2-ARB1)/DLOG(EEND/ESTART)
         B=ARB1-SLOPE*DLOG(ESTART)
         TERM=SLOPE*DLOG(ENER)+B
         SPHI=EXP(TERM)
      ELSEIF(ITY.EQ.9) THEN
         SDEV=EINDX2*CFWSD
         ZZ=(ENER-EINDX1)/SDEV
         SPHI=EXP(-ZZ**2)
      ELSEIF(ITY.EQ.10) THEN
         IF(ENER.LT.0.53)CROSS=0.84*ENER**(-2.81)
         IF(ENER.GE.0.53)CROSS=2.11*ENER**(-2.6)
         OPAC=10.0**(EINDX1-22)*CROSS
         SPHI=0.
         IF(OPAC.LT.100)SPHI=EXP(-OPAC)
         RETURN
       ELSEIF(ITY.EQ.11.OR.ITY.EQ.12) THEN
C        IF(ITY.EQ.11)EDGE=FEEDGE
C        IF(ITY.EQ.12)EDGE=EINDX2
         EDGE=FEEDGE
         ENE=ENER
         IF(ITY.EQ.12)ENE=ENER*FEEDGE/EINDX2
         IF(ENE.LT.EDGE)CROSS=SECTON(ENE,FELABS)
         IF(ENE.GE.EDGE)CROSS=SECTON(ENE,FEKABS)
C        OPAC=CROSS*10.0**(EINDX1-24) <CHANGE ASSUMES ALLEN 6/11>
         OPAC=CROSS*10.0**(EINDX1-4.4-24)
         SPHI=0.
         IF(OPAC.LT.100)SPHI=EXP(-OPAC)
         RETURN
C************* KITAMOTO CHANGE ***********************
C                     REFERENCE BASCO,1980,A&A,87,330
C                     --> BURGESS,1964,MEM.ROY.AST.SOC.,69,1
C  NOTE] HE-LIKE TO H-LIKE NO KUBETSU WA EEKAGEN
       ELSEIF(ITY.EQ.13.OR.ITY.EQ.14) THEN
         SPHI=1.0
         EDGE=EINDX2
         IF(EINDX2.LE.0.02) THEN
           EDGE=FE26ED
           IF(ITY.EQ.14)EDGE=FE25ED
         ENDIF
         ENE=ENER
         IF(ENE.LT.EDGE)RETURN
         CROSS=6.3E-18/26./26.*(EDGE/ENE)**2.84
         IF(ITY.EQ.14)CROSS=CROSS*2
         OPAC=CROSS*10.0**(EINDX1-4.4)
C        OPAC=CROSS*10.0**(EINDX1-4.4)  --- ASSUMES ALLEN
         SPHI=0.
         IF(OPAC.LT.100)SPHI=EXP(-OPAC)
         RETURN
C***********************************  EEKAGEN ********
C     ELSEIF(ITY.EQ.13) THEN
C        SPHI=1.
C        IF(EINDX2.LE.0) EINDX2=FE26ED
C        ENE=ENER*FE26ED/EINDX2
C        IF(ENE.LT.FE26ED)RETURN
C        CROSS=(2./(26.*26.))*(FE26ED/ENE)**2.6666
C        OPAC=CROSS*10.0**(EINDX1-17.)
C        SPHI=0.
C        IF(OPAC.LT.100)SPHI=EXP(-OPAC)
C        RETURN
C     ELSEIF(ITY.EQ.14) THEN
C        SPHI=1.
C        IF(EINDX2.LE.0) EINDX2=FE25ED
C        ENE=ENER*FE25ED/EINDX2
C        IF(ENE.LT.FE25ED)RETURN
C        CROSS=(1.3/(26.*26.))*(FE25ED/ENE)**2.6666
C        OPAC=CROSS*10.0**(EINDX1-17.)
C        SPHI=0.
C        IF(OPAC.LT.100)SPHI=EXP(-OPAC)
C        RETURN
C*************************************************
      ELSEIF(ITY.EQ.15) THEN
         SPHI=1.
         IF(ENER.LE.EINDX1)RETURN
         SPHI=EXP((EINDX1-ENER)/EINDX2)
         RETURN
      ELSEIF(ITY.EQ.16) THEN
ce       ENERA=ALOG10(ENER)
         ENERA=DLOG10(ENER)
         DO 171 I=2,NUMBKG
         IF(ENERA.LT.CBACK(I,2))GO TO 174
171      CONTINUE
         I=NUMBKG
174      I=I-1
         SLOPE=(CBACK(I,1)-CBACK(I+1,1))/(CBACK(I,2)-CBACK(I+1,2))
         TERM2=CBACK(I+1,1)+(ENERA-CBACK(I+1,2))*SLOPE
         SPHI=10.**TERM2
CC
CC  LINE ABSORPTION -- ITY=20: GAUSSIAN OPACITY (DOPPLER)
CC                  -- ITY=21: LORENTZIAN OPACITY (DAMPING)
CC                  -- ITY=22: GAUSSIAN   ABS (TO GET EQUIVALENT WIDTH)
CC                  -- ITY=23: LORENTZIAN ABS (TO GET EQUIVALENT WIDTH)
CC                  -- ITY=24: TRIANGULAR ABS (TO GET EQUIVALENT WIDTH)
CC   DIPOLE= 2.8E-21 KEV*CM**2
CC
      ELSEIF(ITY.EQ.20.OR.ITY.EQ.21.) THEN
         DE=ENER-EINDX2
         IF(ITY.EQ.20)THEN
           DE2=DE*DE/EINDX3/EINDX3
           PROFIL=0.
           IF(DE2.LT.100) PROFIL=EXP(-DE2)/PISQ/EINDX3
         ELSEIF(ITY.EQ.21)THEN
           G=EINDX3/4./PI
           PROFIL=G/(DE**2+G**2)/PI
         ENDIF
         OPAC=DIPOLE*PROFIL*10**(EINDX1-21)
         SPHI=0.
         IF(OPAC.LT.100)SPHI=EXP(-OPAC)
         RETURN
      ELSEIF(ITY.EQ.22)THEN
        DE2=(ENER-EINDX2)/2./EINDX3
        DE2=DE2*DE2
        SPHI=0
        IF(DE2.LT.100) SPHI=EINDX1/1.414213562/PISQ/EINDX3*EXP(-DE2)
        SPHI=1.-SPHI
        IF(SPHI.LT.0) SPHI=0.
        RETURN
      ELSEIF(ITY.EQ.23)THEN
        DE=ENER-EINDX2
        G=EINDX3/4./PI
        SPHI=EINDX1/PI*G/(DE*DE+G*G)
        SPHI=1.-SPHI
        IF(SPHI.LT.0) SPHI=0.
        RETURN
      ELSEIF(ITY.EQ.24) THEN
        IF(EINDX1.GT.1.) EINDX1=1.
        IF(EINDX1.LT.0.) EINDX1=0.
        DE=ABS(ENER-EINDX2)
        SPHI=1.
        IF(DE.LT.EINDX3) SPHI=(1.-EINDX1)/EINDX3*DE+EINDX1
        RETURN
      ELSEIF(ITY.EQ.30) THEN
C******* FIREMAN OPACITY ***84.1.31   S.K.
        CALL FIREMA(ENER,EINDX1,SPHI)
         RETURN
      ELSEIF(ITY.EQ.31) THEN
C******* USER  OPACITY ***84.2. 8   S.K.
        CALL USERAB(ENER,EINDX1,SPHI)
         RETURN
      ELSEIF(ITY.EQ.32) THEN
C******* A*E**(-ALPH)  ***84.2. 8   S.K.
        CALL AEALPH(ENER,EINDX1,EINDX2,EINDX3,SPHI)
         RETURN
      ELSEIF(ITY.EQ.33) THEN
C******* THOMSON SCATTERING **84.2.27   S.K.
        CALL THOMS(ENER,EINDX1,SPHI)
         RETURN
      ELSEIF(ITY.EQ.34) THEN
C******* THOMSON SCATTERING  E DIPENDENT**84.3.12   S.K.
        CALL THOMS2(ENER,EINDX1,SPHI)
         RETURN
      ELSEIF(ITY.EQ.35.OR.ITY.EQ.36) THEN
C******* COMPTONIZED B-B *********  ******84.8.6  S.K./9.4. S.K.
C      EINDX1=KT(B.B) KEV (FOR 35)/ ALPH (FOR36 )
C      EINDX2=TAU     OPTICAL DEPTH
C      EINDX3=KT(ELECTRON) KEV
C
         PROB=1.-EXP(-EINDX2)
         GAMMA=1.0+7.8278E-3*EINDX3+4.5955E-5*EINDX3**2
         NPORI=20
         CALL BCOMP(EINDX1,GAMMA,PROB,ENER,NPORI,FOUT,ITY)
         SPHI=FOUT
      ELSEIF(ITY.EQ.38.OR.ITY.EQ.39) THEN
C******* COMPTONIZED B-B (OCULTATION)*****84.11.2  S.K./9.4. S.K.
C      EINDX1=KT(B.B) KEV (FOR 38)/ ALPH (FOR39 )
C      EINDX2=TAU     OPTICAL DEPTH
C      EINDX3=KT(ELECTRON) KEV
C
         PROB=1.-EXP(-EINDX2)
         GAMMA=1.0+7.8278E-3*EINDX3+4.5955E-5*EINDX3**2
         NPORI=20
         CALL BCOMPO(EINDX1,GAMMA,PROB,ENER,NPORI,FOUT,ITY)
         SPHI=FOUT
      ELSEIF(ITY.EQ.45) THEN
C******* Morrison and McCammon *** 1986/9/18 S.K.
C      EINDX1= column density
C      ENER  = energy
C
         CALL MANDM(ENER,CROSS)
         IF(EINDX1.GT.60.)EINDX1=60.
         OPAC=10.0**EINDX1*CROSS
         SPHI=0.
         IF(OPAC.LT.60)SPHI=EXP(-OPAC)
         RETURN
      ELSEIF(ITY.EQ.46) THEN
C******* Morrison and McCammon/partial absorption *** 1996/6/11 S.K.
C      EINDX1= column density
C      ENER  = energy
C      eindx2: covering factor
C
         CALL MANDM(ENER,CROSS)
         IF(EINDX1.GT.60.)EINDX1=60.
         OPAC=10.0**EINDX1*CROSS
         SPHI=1.-eindx2
         IF(OPAC.LT.60)SPHI=1.-eindx2+eindx2*EXP(-OPAC)
         RETURN
      ENDIF
150   SPHI=SPHI*SCALE
      RETURN
      END
      BLOCK DATA
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON /COEF/ COEFB(7,2),GAUSS(5,2),BREAK(7),CBACK(15,2)
      COMMON /NCOEF/ NUMGU,NUMBKG
      DATA NUMGU /5/
      DATA COEFB /2.,1.9066,1.74,1.59,1.47,1.30,1.,
     *            .7,.74666,0.88,1.08,1.32,2.00,5./
      DATA GAUSS /.236926885056189,.478628670499366,.568888888888888,
     *            .478628670499366,.236926885056189,
     *           -.906179845938664,-.538469310105683,0.0,
     *            .538469310105683, .906179845938664/
      DATA BREAK /1.E38,2.,1.25,.75,.5,.25,.1/
      DATA NUMBKG/7/
      DATA CBACK /.5,-.5,-1.5,-2.5,-3.85,-4.9,-5.65,8*0.,
     *            0.,1.,1.6,2.,2.47,3.,3.47,8*0./
      END
      FUNCTION SECTON(ENER,COEFAT)
      IMPLICIT REAL*8 (A-H,O-Z)
C
C     THE PURPOSE OF SECTON IS TO EVALUATE AN CROSS SECTION FORMULA
C
      DIMENSION COEFAT(4)
ce    ELNB=ALOG(ENER)
      ELNB=DLOG(ENER)
      PSUM=0.
      DO 118 I=1,4
      ELN=ELNB**(I-1)
      PSUM=PSUM+ELN*COEFAT(I)
118   CONTINUE
      SECTON=EXP(PSUM)
      RETURN
      END
C**********************************************************************
C        CHAPLINE & STEVENS (A.P.J. 184 1041(1973))
C     COMPUTON SCATTERING NI YOTE BREMSSTRAHRUNG KARA ZURERU
C        KINGI SHIKI
C          ELECTRON TEMPERATURE  ;T
C          THOMSON SCATTERING LENGTH  ;R
C          PHOTON ENERGY  ;  E
C**********************************************************************
      SUBROUTINE CHASTE(E,T,R,CR)
      IMPLICIT REAL*8 (A-H,O-Z)
C----------------------------------------------------------------------
C     INPUT ARGUMENTS
C        E  :  PHOTON ENERGY (KEV)
C        T  :  ELECTRON TENPERATURE (KEV)
C        R  :  THOMSON SCATTERING LENGTH
C     OUTPUT AREGUMENT
C        CR :  CRECTION FACTER OF FREE FREE RADIATION
C
C        CR=1+B-C+D
C            A=2*T/510
C            ARS=A*R**2
C            EKT=E/T ,EKT2=(E/T)**2 ,EKT3=(E/T)**3
C            P=3*(EKT3-EKT2-EKT)/(2*ARS)
C            Q=3*(2*EKT2-EKT3)/(4*ARS)
C           B=3*EKT3*ARS/8
C           C=P*(EXP(-ARS)-1+ARS)
C           D=Q*(EXP(-ARS)*(ARS+1)-1)
C-----------------------------------------------------------------
             A=2.0*T/510.0
             ARS=A*R**2.0
              IF(ARS.LT.0.000002) THEN
                CR=1.0
                RETURN
              ENDIF
             EKT=E/T
             EKT2=(E/T)**2
             EKT3=(E/T)**3
             P=3.0*(EKT3-EKT2-EKT)/(2.0*ARS)
             Q=3.0*(2.0*EKT2-EKT3)/(4.0*ARS)
            B=3.0*EKT3*ARS/8.0
            C=P*(EXP(-ARS)-1.0+ARS)
            D=Q*(EXP(-ARS)*(ARS+1.0)-1.0)
         CR=1+B-C+D
         RETURN
         END
C**********************************************************
        SUBROUTINE FIREMA(ENER,EINDX1,SPHI)
      IMPLICIT REAL*8 (A-H,O-Z)
C---------------------------------------------------------------------
C        FIREMAN NO CROSS SECTION O TSUKAU  (FIREMAN 1974,AP.J.,187,57
C           GAS NO TOKI
C            ELEMENT  O,NE,MG,SI,S,AR,CA,CR,FE,NI
C
C------------------------------------------------------------------
         DIMENSION EDGE(11),COEF(11),EIND(11)
C----------- K-EDGE ---------------------------------------------
         DATA EDGE/  0.30,  0.53,  0.87,  1.30,  1.84,  2.47,
     &               3.20,  4.04,  5.99,  7.11,  8.33/
C----------- INDEX ----------------------------------------------
         DATA EIND/ -2.82, -2.90, -2.90, -2.85, -2.70, -2.90,
     &              -2.80, -2.81, -2.72, -2.65, -2.76/
C----------- KEISU *10**23 --------------------------------------
         DATA COEF/  9.53, 20.25, 24.34, 26.44, 26.65, 35.62,
     &              34.08, 35.28, 30.59, 52.85, 70.34/
C-----------------------------------------------------------------
          IF(ENER.LT.EDGE(1)) THEN
           CROSS=2.11*ENER**(-2.6)
           GO TO 3000
          ENDIF
 1001     DO 1000 I1=1,10
          I1A=I1
         IF(ENER.GE.EDGE(I1).AND.ENER.LT.EDGE(I1+1))GO TO 2000
 1000     CONTINUE
          I1A=11
 2000     CONTINUE
            CROSS=COEF(I1A)*ENER**(EIND(I1A))
 3000     OPAC=10.0**(EINDX1-23)*CROSS
           SPHI=0.0
         IF(OPAC.LT.100)SPHI=EXP(-OPAC)
         RETURN
         END
C**********************************************************
        SUBROUTINE USERAB(ENER,EINDX1,SPHI)
      IMPLICIT REAL*8 (A-H,O-Z)
C---------------------------------------------------------------------
C        USER  NO CROSS SECTION DE KYUSHU KEISU O TSUKURU
C
C       $  ** 84/2/9  .... FOR  PROF. MIYAMOTO
C       X  ** 84/2/14 .... FOR  PARTIAL IONIZED GASS
C                           C(V),N(VI),O(VII),SI(IX),S(IX),FE(XII)
C------------------------------------------------------------------
         DIMENSION EDGE(11),COEF(11),EIND(11)
C----------- K-EDGE ---------------------------------------------
         DATA EDGE/  1.50,  2.00,  3.00,  4.00,  5.00,  7.10,
     &               99.0,  99.0,  99.0, 99.00, 99.00/
C----------- INDEX ----------------------------------------------
         DATA EIND/  1.68, -1.33, -2.51, -3.38,-14.47, -2.65,
     &              -0.00, -0.00, -0.00, -0.00, -0.00/
C----------- KEISU *10**23 --------------------------------------
         DATA COEF/ 0.375,  3.02, 11.03, 36.75,2.09E9, 36.05,
     &              99.99, 99.99, 99.99, 99.99, 99.99/
C-----------------------------------------------------------------
          IF(ENER.LT.EDGE(1)) THEN
           CROSS=0.30*ENER**(2.23)
           GO TO 3000
          ENDIF
 1001     DO 1000 I1=1,10
          I1A=I1
         IF(ENER.GE.EDGE(I1).AND.ENER.LT.EDGE(I1+1))GO TO 2000
 1000     CONTINUE
          I1A=11
 2000     CONTINUE
            CROSS=COEF(I1A)*ENER**(EIND(I1A))
 3000     OPAC=10.0**(EINDX1-23)*CROSS
           SPHI=0.0
         IF(OPAC.LT.100)SPHI=EXP(-OPAC)
         RETURN
         END
C****************************************************
       SUBROUTINE AEALPH(ENER,EINDX1,EINDX2,EINDX3,SPHI)
      IMPLICIT REAL*8 (A-H,O-Z)
C---------------------------------------------------
C        OPECITY O A*E**(-ALPH) TO SURU
C               EINDX1...... NH
C               EINDX2...... A
C               EINDX3...... ALPH
C----------------------------------------------------
          CROSS=EINDX2*ENER**(-EINDX3)
         OPAC=10.0**(EINDX1-22)*CROSS
         SPHI=0.
         IF(OPAC.LT.100)SPHI=EXP(-OPAC)
         RETURN
         END
C*****************************************************
        SUBROUTINE THOMS(ENER,EINDX1,SPHI)
      IMPLICIT REAL*8 (A-H,O-Z)
C-----------------------------------------------------
C        NH ATARI NH NO ELECTRON
C             ELECTRON 1 KO ATARI .... 6.65*10**(-25) CM**2
C             ELECTRON  10**22 ATARI...6.65*10**(-3)  CM**2
C
C------------------------------------------------------
        DATA THO22 /0.00665/
         OPAC=10**(EINDX1-22)*THO22
C        OPAC=10**(EINDX1-22)*THO22*1.2
         SPHI=0.
         IF(OPAC.LT.100)SPHI=EXP(-OPAC)
         RETURN
         END
C*****************************************************
        SUBROUTINE THOMS2(ENER,EINDX1,SPHI)
      IMPLICIT REAL*8 (A-H,O-Z)
C-----------------------------------------------------
C        NH ATARI NH NO ELECTRON
C         PHAI0=  ELECTRON 1 KO ATARI .... 6.65*10**(-25) CM**2
C              =  ELECTRON  10**22 ATARI...6.65*10**(-3)  CM**2
C        W.HEITLER P.226
C            PHAI=PHAI0*(1-2.0*GAMMA+5.2*GAMMA**2)
C                 GAMMA=ENER(KEV)/512(KEV)
C------------------------------------------------------
        DATA THO22 /0.00665/
        DATA ELEC0 /512.0/
         GAMMA=ENER/ELEC0
         OPAC=10**(EINDX1-22)*THO22*(1.0-2.0*GAMMA+5.2*GAMMA**2)
         SPHI=0.
         IF(OPAC.LT.100)SPHI=EXP(-OPAC)
         RETURN
         END
C***************************************************************
       SUBROUTINE BCOMP(T,G,P,E,N,FOUT,ITY)
       IMPLICIT REAL*8 (A-H,O-Z)
C--------------------------------------------------------------
C      COMPTONIZED B-B NO SPECTRUM O MOTOMERU
C       INPUT AREGUMENTS
C          T  : B-B TEMP.(KEV)FOR 35/ ALPH FOR 36
C          G  : (1+4*KTE/M/C**2+12*(KTE/M/C**2)**2) KTE: ELECTRON TE-
C               1.0+7.8278E-3*KTE+4.5955E-5*KTE**2    MPERATURE(KEV)
C          P  : PROBABIRITY OF SCATTERING (1-EXP(-TAU) TAU:OPTICAL
C                                                          DEPTH)
C          E  : HOSHI ENERGY (KEV)
C          N  : PORINOMIAL NO KAISUU
C          ITY: INPUT SPECTRUM TYPE  35 : BB/ 36 : PL
C       OUTPUT AREGUMENTS
C          FOUT : ENERGY E NO PHOTON SUU
C
C        TAU<1.0 DE NAITO TSUKAENAI
C                                        N.SHIBAZAKI & K.MITSUDA
C                              SUBROUTINE: KITAMOTO
C-------------------------------------------------------------
        BB(T,EGN)=EGN*EGN/(EXP(EGN/T)-1.)
        PL(T,EGN)=EGN**(-T)
C--------------- TAU CHECK -----------------
         IF(P.GT.1.0)GO TO 999
C--------------- ITY CHECK -----------------
         IF(ITY.NE.35.AND.ITY.NE.36)GO TO 998
C--------------------------------------------
   1    CONTINUE
        FOUT=0.
 1001   DO 1000 I1=0,N
         FI1=I1
         EGN=E/G**FI1
         IF(EGN.LT.0.000001)GO TO 9000
         IF(ITY.EQ.35)FIN=BB(T,EGN)
         IF(ITY.EQ.36)FIN=PL(T,EGN)
         FOUT=FOUT+FIN*(P/G)**FI1*(1.-P)
 1000   CONTINUE
 9000   RETURN
C--------------- MADA DEKITENAI -----------------
 998      PRINT *,'??????? ITY = ',ITY,'  IN BCOMP ?????'
           STOP
C--------------- TAU G.T. 1 ---------------------
 999      PRINT *,'??????? TAU G.T.1 IN BCOMP ??????????'
          PRINT *,'??????? SET TAU.EQ.1 BCOMP ??????????'
C         P=1.0
          GO TO 1
C------------------------------------------------
        END
C***************************************************************
       SUBROUTINE BCOMPO(T,G,P,E,N,FOUT,ITY)
       IMPLICIT REAL*8 (A-H,O-Z)
C--------------------------------------------------------------
C      COMPTONIZED B-B NO SPECTRUM O MOTOMERU
C       INPUT AREGUMENTS
C          T  : B-B TEMP.(KEV)FOR 35/ ALPH FOR 36
C          G  : (1+4*KTE/M/C**2+12*(KTE/M/C**2)**2) KTE: ELECTRON TE-
C               1.0+7.8278E-3*KTE+4.5955E-5*KTE**2    MPERATURE(KEV)
C          P  : PROBABIRITY OF SCATTERING (1-EXP(-TAU) TAU:OPTICAL
C                                                          DEPTH)
C          E  : HOSHI ENERGY (KEV)
C          N  : PORINOMIAL NO KAISUU
C          ITY: INPUT SPECTRUM TYPE  35 : BB/ 36 : PL
C       OUTPUT AREGUMENTS
C          FOUT : ENERGY E NO PHOTON SUU
C
C        TAU<1.0 DE NAITO TSUKAENAI
C                                        N.SHIBAZAKI & K.MITSUDA
C                              SUBROUTINE: KITAMOTO
C-------------------------------------------------------------
        BB(T,EGN)=EGN*EGN/(EXP(EGN/T)-1.)
        PL(T,EGN)=EGN**(-T)
C--------------- TAU CHECK -----------------
         IF(P.GT.1.0)GO TO 999
C--------------- ITY CHECK -----------------
         IF(ITY.NE.38.AND.ITY.NE.39)GO TO 998
C--------------------------------------------
   1    CONTINUE
        FOUT=0.
 1001   DO 1000 I1=1,N
         FI1=I1
         EGN=E/G**FI1
         IF(EGN.LT.0.000001)GO TO 9000
         IF(ITY.EQ.38)FIN=BB(T,EGN)
         IF(ITY.EQ.39)FIN=PL(T,EGN)
         FOUT=FOUT+FIN*(P/G)**FI1*(1.-P)
 1000   CONTINUE
 9000   RETURN
C--------------- MADA DEKITENAI -----------------
 998      PRINT *,'??????? ITY = ',ITY,'  IN BCOMPO?????'
           STOP
C--------------- TAU G.T. 1 ---------------------
 999      PRINT *,'??????? TAU G.T.1 IN BCOMPO??????????'
          PRINT *,'??????? SET TAU.EQ.1 BCOMPO??????????'
C         P=1.0
          GO TO 1
C------------------------------------------------
        END
C***************************************************************
       SUBROUTINE MANDM(E,CROSS)
       IMPLICIT REAL*8 (A-H,O-Z)
C--------------------------------------------------------------
C      Inter stellar absorption cross section by Morrison and McCammon
C                          (1983) Astrophys J. 270, 119
C       Input Areguments
C          E    : photon energy (keV)
C       Output Areguments
C          CROSS: cross ection in cm**2
C   Note]   This formula valids in the energy range 0.030-10.000 keV.
C
C    Cross section per hydrogen atom
C      =(C0+C1*E+C2*E**2)*E**(-3)*10**(-24)cm**2 (E in keV)
C                              Subroutine: S.Kitamoto
C-------------------------------------------------------------
      DIMENSION C0(14), C1(14), C2(14),ECHAN(15)
      DATA C0 /  17.3,  34.6,  78.1,  71.4,  95.5, 308.9, 120.6,
     &          141.3, 202.7, 342.7, 352.2, 433.9, 629.0, 701.2/
      DATA C1 / 608.1, 267.9,  18.8,  66.8, 145.8,-380.6, 169.3,
     &          146.8, 104.7,  18.7,  18.7,  -2.4,  30.9,  25.2/
      DATA C2 /-2150.,-476.1,   4.3, -51.4, -61.1, 294.0, -47.7,
     &          -31.5, -17.0,   0.0,   0.0,  0.75,   0.0,   0.0/
      DATA ECHAN / 0.030, 0.100, 0.284, 0.400, 0.532,
     &             0.707, 0.867, 1.303, 1.840, 2.471,
     &             3.210, 4.038, 7.111, 8.331,10.000/
 1001 DO 1000 I1=1,14
      IF(E.LT.ECHAN(I1+1)) THEN
       ICHAN=I1
       GO TO 1100
      ENDIF
 1000 CONTINUE
      ICHAN=14
 1100 CROSS = (C0(ICHAN)+C1(ICHAN)*E+C2(ICHAN)*E**(2.))
     &        *E**(-3.)*10.**(-24.)
      RETURN
      END
