C             DEBUG                                                   
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*        RECOIL  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                             00270000
            POWER=1.-EINDX1                                             00280000
            PLCNST=1./POWER                                             00290000
C-------------- KITAMOTO CHANGE --------                                00300001
             IF(POWER.GT.250.)POWER=252.                                00310001
            SUM1=E1**(POWER)                                            00320000
            SUM2=E2**(POWER)                                            00330000
            PHOTON=PLCNST*(SUM2-SUM1)                                   00340000
         ELSE                                                           00350000
            PHOTON=DLOG(E2/E1)                                          00360000
         ENDIF                                                          00370000
C                                                                       00380000
C     DO LIANG THERMAL CYCLTRON                                         00390000
C                                                                       00400000
      ELSEIF(ITY.EQ.2) THEN                                             00410000
         DB=(E2-E1)/NUMINT                                              00420000
         DDB=DB/2.                                                      00430000
         SUM=0.                                                         00440000
         DO 211 I=1,NUMINT                                              00450000
         APBB=(2*I-1)*DDB+E1                                            00460000
         DO 211 II=1,NUMGU                                              00470000
         E=DDB*GAUSS(II,2)+APBB                                         00480000
         AE=4.5*E*EINDX1                                                00490000
         TERM=EXP(-AE**(1./3.))                                         00500000
         SUM=SUM+GAUSS(II,1)*TERM                                       00510000
211      CONTINUE                                                       00520000
         PHOTON=SUM*DDB                                                 00530000
C                                                                       00540000
C     DO THERMAL BREMSSTRAHLUNG (CLASSICAL)                             00550000
C                                                                       00560000
      ELSEIF(ITY.EQ.3) THEN                                             00570000
         DB=(E2-E1)/NUMINT                                              00580000
         DDB=DB/2.                                                      00590000
         SUM=0.                                                         00600000
         DO 32 I=1,NUMINT                                               00610000
         APBB=(2*I-1)*DDB+E1                                            00620000
         DO 32 II=1,NUMGU                                               00630000
         E=DDB*GAUSS(II,2)+APBB                                         00640000
         AE=E/EINDX1                                                    00650000
C---------  CHAPLIN STEVENS SPECTRUM  .... S.K. --------                00660000
           IF(EINDX2.GT.0.0001)THEN                                     00670000
            CALL CHASTE(E,EINDX1,EINDX2,CR)                             00680000
            TERM=TERM*CR                                                00690000
           ENDIF                                                        00700000
C--------------------------------------------------------               00710000
         TERM=EXP(-AE)/E                                                00720000
         SUM=SUM+GAUSS(II,1)*TERM                                       00730000
32       CONTINUE                                                       00740000
         PHOTON=SUM*DDB                                                 00750000
C                                                                       00760000
C     DO BLACK BODY LAW                                                 00770000
C                                                                       00780000
      ELSEIF(ITY.EQ.4) THEN                                             00790000
C        BBCNST=EINDX1**(3.)                                            00800000
C        ARG1=E1/EINDX1                                                 00810000
C        C11=EXP(-ARG1)                                                 00820000
C        EXPQ1=1.                                                       00830000
C        ARG=E2/EINDX1                                                  00840000
C        C1=EXP(-ARG)                                                   00850000
C        EXPQ=1.                                                        00860000
C        SUM1=0.                                                        00870000
C        SUM2=0.                                                        00880000
C        DO 45 N=1,10                                                   00890000
C        Q1=ARG1*N                                                      00900000
C        Q21=Q1*Q1                                                      00910000
C        POLY1=Q21+2.*Q1+2.                                             00920000
C        EXPQ1=C11*EXPQ1                                                00930000
C        SUM1=SUM1+EXPQ1*POLY1/(N**3)                                   00940000
C        Q=ARG*N                                                        00950000
C        Q2=Q*Q                                                         00960000
C        POLY=Q2+2.*Q+2.                                                00970000
C        EXPQ=C1*EXPQ                                                   00980000
C        SUM2=SUM2+EXPQ*POLY/(N**3)                                     00990000
C45       CONTINUE                                                      01000000
C         PHOTON=BBCNST*(SUM1-SUM2)                                     01010000
         DB=(E2-E1)/NUMINT                                              01020000
         DDB=DB/2.                                                      01030000
         SUM=0.                                                         01040000
         DO 46 I=1,NUMINT                                               01050000
         APBB=(2*I-1)*DDB+E1                                            01060000
         DO 46 II=1,NUMGU                                               01070000
         E=DDB*GAUSS(II,2)+APBB                                         01080000
         U=E/EINDX1                                                     01090000
C----------- KITAMOTO CHANGE ------------                               01100000
C        IF(U.GT.500.) U=500 --- X.LT.174.673 FOR EXP(X)                01110000
         IF(U.GT.170.) GO TO 46                                         01120000
         TERM=E*E/(EXP(U)-1.)                                           01130000
         SUM=SUM+GAUSS(II,1)*TERM                                       01140000
46       CONTINUE                                                       01150000
         PHOTON=SUM*DDB                                                 01160000
C                                                                       01170000
C     DO THERM BREM, MATTESON OR COMBINED MATTESON-JACOBSON             01180000
C                                                                       01190000
      ELSEIF(ITY.EQ.5.OR.ITY.EQ.6) THEN                                 01200000
         UNIT=EINDX1/100.                                               01210000
         PW1=-.333                                                      01220000
         IF((1./EINDX1).GT..01)PW1=-.190                                01230000
         PW2=-.31*UNIT**PW1                                             01240000
         PW3=.415                                                       01250000
         IF((1./EINDX1).GT..01)PW3=.125                                 01260000
         COEFTB=1.04*UNIT**PW3                                          01270000
         IF(ITY.EQ.6) THEN                                              01280000
            DO 69 I=1,7                                                 01290000
            COEFBU(I,1)=1.+COEFB(I,1)*EINDX1/511.                       01300000
            COEFBU(I,2)=COEFB(I,2)/511.                                 01310000
69          CONTINUE                                                    01320000
         ENDIF                                                          01330000
         DB=(E2-E1)/NUMINT                                              01340000
         DDB=DB/2.                                                      01350000
         SUM=0.                                                         01360000
         DO 51 I=1,NUMINT                                               01370000
         APBB=(2*I-1)*DDB+E1                                            01380000
         DO 51 II=1,NUMGU                                               01390000
         E=DDB*GAUSS(II,2)+APBB                                         01400000
         AE=E/EINDX1                                                    01410000
C------------ KITAMOTO ----------------                                 01420001
           IF(AE.GT.170.)AE=170.                                        01430001
         TERM=EXP(-AE)*AE**PW2/E                                        01440000
         IF(ITY.EQ.6)THEN                                               01450000
            IBST=7                                                      01460000
            DO 60 WHILE(AE.GT.BREAK(IBST))                              01470000
60          IBST=IBST-1                                                 01480000
            TERM=TERM*(COEFBU(IBST,1)+COEFBU(IBST,2)*E)                 01490000
         ENDIF                                                          01500000
C---------  CHAPLIN STEVENS SPECTRUM  .... S.K. --------                01510000
           IF(EINDX2.GT.0.0001)THEN                                     01520000
            CALL CHASTE(E,EINDX1,EINDX2,CR)                             01530000
            TERM=TERM*CR                                                01540000
           ENDIF                                                        01550000
C--------------------------------------------------------               01560000
         SUM=SUM+GAUSS(II,1)*TERM                                       01570000
51       CONTINUE                                                       01580000
         PHOTON=SUM*COEFTB*DDB                                          01590000
C                                                                       01600000
C     DO EXPONENTIAL                                                    01610000
C                                                                       01620000
      ELSEIF(ITY.EQ.7) THEN                                             01630000
         SUM1=EXP(-E1/EINDX1)*EINDX1                                    01640000
         SUM2=EXP(-E2/EINDX1)*EINDX1                                    01650000
         PHOTON=SUM1-SUM2                                               01660000
C                                                                       01670000
C     DO SEGMENT OF A POWER LAW DEFINED BY TWO POINTS                   01680000
C                                                                       01690000
      ELSEIF(ITY.EQ.8) THEN                                             01700000
         ARB1=DLOG(EINDX1)                                              01710000
         ARB2=DLOG(EINDX2)                                              01720000
         SLOPE=(ARB2-ARB1)/DLOG(EEND/ESTART)                            01730000
         BB=ARB1-SLOPE*DLOG(ESTART)                                     01740000
         SLOPEP=SLOPE+1.                                                01750000
         IF(ABS(SLOPEP).LE.(.01)) THEN                                  01760000
            PHOTON=EXP(BB)*DLOG(E2/E1)                                  01770000
         ELSE                                                           01780000
            SUM=(E2**SLOPEP-E1**SLOPEP)/SLOPEP                          01790000
            PHOTON=SUM*EXP(BB)                                          01800000
         ENDIF                                                          01810000
C                                                                       01820000
C     DO GAUSSIAN LINE, NOTE THAT LOADFR ASSUME THAT GAUSSIAN IS        01830000
C     FORMULA NUM 9                                                     01840000
C                                                                       01850000
      ELSEIF(ITY.EQ.9) THEN                                             01860000
         DB=(E2-E1)/NUMINT                                              01870000
         DDB=DB/2.                                                      01880000
         SUM=0.                                                         01890000
         W=CFWSD*EINDX2                                                 01900000
         DO 102 I=1,NUMINT                                              01910000
         APBB=(2*I-1)*DDB+E1                                            01920000
         DO 102 II=1,NUMGU                                              01930000
         E=DDB*GAUSS(II,2)+APBB                                         01940000
         ZZ=(E-EINDX1)/W                                                01950000
         TERM=EXP(-ZZ**2)                                               01960000
         SUM=SUM+GAUSS(II,1)*TERM                                       01970000
102      CONTINUE                                                       01980000
         PHOTON=SUM*DDB/(W*1.77245385)                                  01990000
C                                                                       02000000
C     SUMMING OF VARIOUS ABSORPTIONS IS NOT ALLOWED, IT ONLY IS USED    02010000
C     FOR MULTIPLICATION (IN PHIV)                                      02020000
C                                                                       02030000
      ELSEIF(ITY.GE.10.AND.ITY.LE.14) THEN                              02040000
         PHOTON=0                                                       02050000
C******* FOR FIREMAN ABSORP TION 84.1.31  S.K. *******                  02060000
      ELSEIF(ITY.EQ.30) THEN                                            02070000
         PHOTON=0                                                       02080000
C******* FOR USER ABSORPTION     84.1.31  S.K. *******                  02090004
      ELSEIF(ITY.EQ.31) THEN                                            02100004
         PHOTON=0                                                       02110004
C******* A*E**(-ALPH) ABSORPTION 84.1.31  S.K. *******                  02120004
      ELSEIF(ITY.EQ.32) THEN                                            02130004
         PHOTON=0                                                       02140004
C******* A*E**(-ALPH) ABSORPTION 84.1.31  S.K. *******                  02150004
      ELSEIF(ITY.EQ.33) THEN                                            02160004
         PHOTON=0                                                       02170004
C                                                                       02180000
C     DO EXP CUTOFF                                                     02190000
C                                                                       02200000
      ELSEIF(ITY.EQ.15) THEN                                            02210000
         DB=(E2-E1)/NUMINT                                              02220000
         DDB=DB/2.                                                      02230000
         SUM=0.                                                         02240000
         DO 162 I=1,NUMINT                                              02250000
         APBB=(2*I-1)*DDB+E1                                            02260000
         DO 162 II=1,NUMGU                                              02270000
         E=DDB*GAUSS(II,2)+APBB                                         02280000
         TERM=EXP((EINDX1-E)/EINDX2)                                    02290000
         SUM=SUM+GAUSS(II,1)*TERM                                       02300000
162      CONTINUE                                                       02310000
         PHOTON=SUM*DDB                                                 02320000
C                                                                       02330000
C     DO DIFFUSE BACKGROUND, NOTE THAT FORSC SHOULD TO THE STER OF      02340000
C     THE SYSTEM                                                        02350000
C                                                                       02360000
      ELSEIF(ITY.EQ.16) THEN                                            02370000
         DB=(E2-E1)/NUMINT                                              02380000
         DDB=DB/2.                                                      02390000
         SUM=0.                                                         02400000
         DO 172 I=1,NUMINT                                              02410000
         APBB=(2*I-1)*DDB+E1                                            02420000
         DO 172 II=1,NUMGU                                              02430000
         E=DDB*GAUSS(II,2)+APBB                                         02440000
         ENERA=DLOG10(E)                                                02450000
         DO 171 K=2,NUMBKG                                              02460000
         IF(ENERA.LT.CBACK(K,2))GO TO 174                               02470000
171      CONTINUE                                                       02480000
         K=NUMBKG                                                       02490000
174      K=K-1                                                          02500000
         SLOPE=(CBACK(K,1)-CBACK(K+1,1))/(CBACK(K,2)-CBACK(K+1,2))      02510000
         TERM2=CBACK(K+1,1)+(ENERA-CBACK(K+1,2))*SLOPE                  02520000
         TERM=10.**TERM2                                                02530000
         SUM=SUM+GAUSS(II,1)*TERM                                       02540000
172      CONTINUE                                                       02550000
         PHOTON=SUM*DDB                                                 02560000
      ELSEIF(ITY.EQ.18) THEN                                            02570026
C******* DISK MODEL **************  ******85.2.15   S.K.(COPY FROM K.M) 02570126
C        EINDX1 : KT                                                   
         DB=(E2-E1)/NUMINT                                             
         DDB=DB/2.                                                     
         SUM=0.                                                        
         DO 47 I=1,NUMINT                                              
         APBB=(2*I-1)*DDB+E1                                           
         DO 47 II=1,NUMGU                                              
         E=DDB*GAUSS(II,2)+APBB   
*         print *,"eindx1=",eindx1," e=",e
         CALL DSKSPC(EINDX1,E,FLUX)                                    
*         print *,"flux=",flux
         TERM=FLUX                                                     
         SUM=SUM+GAUSS(II,1)*TERM                                      
47       CONTINUE                                                      
         PHOTON=SUM*DDB                                                
      ELSEIF(ITY.EQ.35.OR.ITY.EQ.36) THEN                               02571726
C******* COMPTONIZED B-B *********  ******84.8.6    S.K./9.4 S.K.       02580004
C      EINDX1=KT(B.B) KEV(FOR 35)/ ALPH(P.L.) (FOR.36)                  02590004
C      EINDX2=TAU     OPTICAL DEPTH                                     02600000
C      EINDX3=KT(ELECTRON) KEV                                          02610000
C                                                                       02620000
         PROB=1.-EXP(-EINDX2)                                           02630000
         GAMMA=1.0+7.8278E-3*EINDX3+4.5955E-5*EINDX3**2                 02640000
         NPORI=20                                                       02650004
         DB=E2-E1                                                       02660003
         E=DB/2.+E1                                                     02670003
         CALL BCOMP(EINDX1,GAMMA,PROB,E,NPORI,FOUT,ITY)                 02680003
         PHOTON=FOUT*DB                                                 02690003
*      ELSEIF(ITY.EQ.37) THEN                                            02700004
C******* COMPTONIZED P-L *BY MITSUDA******84.8.6    S.K./9.4 S.K.       02710002
C     COMMON /PARMS/TS,TE,ODEPTH,EF,X                                   02720004
C      EINDX1= ALPH (PHOTON INDEX)                                      02730004
C      EINDX2=TAU     OPTICAL DEPTH                                     02740004
C      EINDX3=KT(ELECTRON) KEV                                          02750004
C                                                                       02760004
*         TS=EINDX1                                                      02770004
*         ODEPTH=EINDX2                                                  02780004
*         TE=EINDX3                                                      02790004
*         EF=(E2+E1)/2.                                                  02800004
*         CALL SETK3                                                     02810004
*         CALL CMPSUN(FOUT)                                              02820004
*         TERM=FOUT                                                      02830004
*         PHOTON=TERM*(E2-E1)                                            02840004
      ELSEIF(ITY.EQ.38.OR.ITY.EQ.39) THEN                               02850002
C******* COMPTONIZED B-B (OCULTATION)*****84.11.2    S.K./9.4 S.K.      02860002
C      EINDX1=KT(B.B) KEV(FOR 38)/ ALPH(P.L.) (FOR.39)                  02870002
C      EINDX2=TAU     OPTICAL DEPTH                                     02880002
C      EINDX3=KT(ELECTRON) KEV                                          02890002
C                                                                       02900002
         PROB=1.-EXP(-EINDX2)                                           02910002
         GAMMA=1.0+7.8278E-3*EINDX3+4.5955E-5*EINDX3**2                 02920002
         NPORI=20                                                       02930002
         DB=E2-E1                                                       02940003
         E=DB/2.+E1                                                     02950003
         CALL BCOMPO(EINDX1,GAMMA,PROB,E,NPORI,FOUT,ITY)                02960002
         PHOTON=FOUT*DB                                                 02970003
      ELSEIF(ITY.EQ.40.OR.ITY.EQ.41.OR.ITY.EQ.42) THEN                  02980011
C******* COMPTONIZED B-B (OCULTATION)*****84.11.5    S.K.               02990004
C      EINDX1=KT(B.B) KEV(FOR 40,42)/ ALPH(P.L.) (FOR.41)               03000011
C      EINDX2=TAU     OPTICAL DEPTH                                     03010004
C      EINDX3=KT(ELECTRON) KEV                                          03020004
C       ITY.EQ.42 : OCULTATION                                          03030011
C                                                                       03040004
         PROB=1.-EXP(-EINDX2)                                           03050004
         F1KTMC=1.+2*EINDX3/511.002                                     03060004
         NPORI=20                                                       03070004
         DB=E2-E1                                                       03080004
         E=DB/2.+E1                                                     03090004
         CALL BCOMPR(EINDX1,F1KTMC,PROB,E,NPORI,FOUT,ITY)               03100004
         PHOTON=FOUT*DB                                                 03110004
      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

      ENDIF                                                             03120004
C                                                                       03130000
C     DO RETURN                                                         03140000
C                                                                       03150000
      PHOTON=PHOTON*SCALE                                               03160000
      RETURN                                                            03170000
      END                                                               03180000
      SUBROUTINE PHIV(EINDX1,EINDX2,ITY,SCALE,ENER,ESTART,EEND,         03190000
     &  EINDX3,EINDX4,SPHI)                                             03200000
      IMPLICIT REAL*8 (A-H,O-Z)                                         03210004
C                                                                       03220000
C     THE PURPOSE OF PHIV IS TO GIVE THE VALUE OF THE PHOTON SPECTRUM   03230000
C     AT ENER                                                           03240000
C                                                                       03250000
      COMMON /PARMS/TS,TE,ODEPTH,EF,X                                   03260004
C      COMMON /COEF/ NUMGU,COEFB(7,2),GAUSS(5,2),BREAK(7),NUMBKG,        03270000
C     *               CBACK(15,2)                                        03280000
      COMMON /COEF/ COEFB(7,2),GAUSS(5,2),BREAK(7), CBACK(15,2)
      COMMON /NCOEF/ NUMGU,NUMBKG

      DATA NUMINT/1/,CFWSD/0.60056D0/                                         03290000
      DIMENSION FELABS(4),FEKABS(4),COEFBU(7,2)                         03300000
      DATA FELABS/13.6696,-2.39195,-.137648,0./                         03310000
      DATA FEKABS/14.3456,-1.23491,-.418785,.0321662/                   03320000
      DATA FEEDGE/7.111/                                                03330000
      DATA FE26ED,FE25ED/9.2,8.7/                                       03340000
C      DATA CFWSD/0.6D0/                                                   03350000
      DATA DIPOLE,PI,PISQ/2.8,3.141592653,1.77245385/                   03360000
      IF(ITY.EQ.1) THEN                                                 03370000
         SPHI=ENER**(-EINDX1)                                           03380000
      ELSEIF(ITY.EQ.2) THEN                                             03390000
         AE=4.5*ENER*EINDX1                                             03400000
         SPHI=EXP(-AE**.333333333)                                      03410000
      ELSEIF(ITY.EQ.3) THEN                                             03420000
         AE=ENER/EINDX1                                                 03430000
         IF(AE.GT.500)AE=500                                            03440000
C----------- CHAPLINE STEVENS SPECTRUM..... S.K. --------               03450000
           CR=1.0000                                                    03460000
           IF(EINDX2.GT.0.0001)CALL CHASTE(ENER,EINDX1,EINDX2,CR)       03470000
           SPHI=EXP(-AE)/ENER *CR                                       03480000
C--------------------------------------------------------               03490000
C        SPHI=EXP(-AE)/ENER                                             03500000
      ELSEIF(ITY.EQ.4) THEN                                             03510000
         U=ENER/EINDX1                                                  03520000
         IF(U.GT.500)U=500                                              03530000
         SPHI=ENER*ENER/(EXP(U)-1.)                                     03540000
      ELSEIF(ITY.EQ.5.OR.ITY.EQ.6) THEN                                 03550000
         UNIT=.01*EINDX1                                                03560000
         PW1=-.3333                                                     03570000
         IF((1./EINDX1).GT..01)PW1=-.190                                03580000
         PW2=-.31*UNIT**PW1                                             03590000
         PW3=.415                                                       03600000
         IF((1./EINDX1).GT..01)PW3=.125                                 03610000
         AE=ENER/EINDX1                                                 03620000
         IF(AE.GT.500)AE=500                                            03630000
         COEFTB=1.04*UNIT**PW3                                          03640000
         SPHI=COEFTB*EXP(-AE)*AE**PW2/ENER                              03650000
         IF(ITY.EQ.6) THEN                                              03660000
            DO 69 I=1,7                                                 03670000
            COEFBU(I,1)=1.+COEFB(I,1)*EINDX1/511.                       03680000
            COEFBU(I,2)=COEFB(I,2)/511.                                 03690000
69          CONTINUE                                                    03700000
            IBST=7                                                      03710000
            DO 61 WHILE(AE.GT.BREAK(IBST))                              03720000
61          IBST=IBST-1                                                 03730000
            SPHI=SPHI*(COEFBU(IBST,1)+COEFBU(IBST,2)*ENER)              03740000
         ENDIF                                                          03750000
C----------- CHAPLINE STEVENS SPECTRUM..... S.K. --------               03760000
           CR=1.0000                                                    03770000
           IF(EINDX2.GT.0.0001)CALL CHASTE(ENER,EINDX1,EINDX2,CR)       03780000
           SPHI=SPHI *CR                                                03790000
C--------------------------------------------------------               03800000
      ELSEIF(ITY.EQ.7) THEN                                             03810000
         AE=ENER/EINDX1                                                 03820000
         IF(AE.GT.500)AE=500                                            03830000
         SPHI=EXP(-AE)                                                  03840000
      ELSEIF(ITY.EQ.8) THEN                                             03850000
         ARB1=DLOG(EINDX1)                                              03860000
         ARB2=DLOG(EINDX2)                                              03870000
         SLOPE=(ARB2-ARB1)/DLOG(EEND/ESTART)                            03880000
         B=ARB1-SLOPE*DLOG(ESTART)                                      03890000
         TERM=SLOPE*DLOG(ENER)+B                                        03900000
         SPHI=EXP(TERM)                                                 03910000
      ELSEIF(ITY.EQ.9) THEN                                             03920000
         SDEV=EINDX2*CFWSD                                              03930000
         ZZ=(ENER-EINDX1)/SDEV                                          03940000
         SPHI=EXP(-ZZ**2)                                               03950000
      ELSEIF(ITY.EQ.10) THEN                                            03960000
         IF(ENER.LT.0.53)CROSS=0.84*ENER**(-2.81)                       03970000
         IF(ENER.GE.0.53)CROSS=2.11*ENER**(-2.6)                        03980000
C@@      IF(EINDX1.GT.70)EINDX1=70                                      03990021
         OPAC=10.0**(EINDX1-22)*CROSS                                   04000000
         SPHI=0.                                                        04010000
         IF(OPAC.LT.100)SPHI=EXP(-OPAC)                                 04020000
         RETURN                                                         04030000
       ELSEIF(ITY.EQ.11.OR.ITY.EQ.12) THEN                              04040000
C        IF(ITY.EQ.11)EDGE=FEEDGE                                       04050000
C        IF(ITY.EQ.12)EDGE=EINDX2                                       04060000
         EDGE=FEEDGE                                                    04070000
         ENE=ENER                                                       04080000
         IF(ITY.EQ.12)ENE=ENER*FEEDGE/EINDX2                            04090000
         IF(ENE.LT.EDGE)CROSS=SECTON(ENE,FELABS)                        04100000
         IF(ENE.GE.EDGE)CROSS=SECTON(ENE,FEKABS)                        04110000
C@@      IF(EINDX1.GT.70)EINDX1=70                                      04120021
C        OPAC=CROSS*10.0**(EINDX1-24) <CHANGE ASSUMES ALLEN 6/11>       04130000
         OPAC=CROSS*10.0**(EINDX1-4.4-24)                               04140000
         SPHI=0.                                                        04150000
         IF(OPAC.LT.100)SPHI=EXP(-OPAC)                                 04160000
         RETURN                                                         04170000
C************* KITAMOTO CHANGE ***********************                  04180004
C                     REFERENCE BASCO,1980,A&A,87,330                   04190004
C                     --> BURGESS,1964,MEM.ROY.AST.SOC.,69,1            04200004
C  NOTE] HE-LIKE TO H-LIKE NO KUBETSU WA EEKAGEN                        04210004
       ELSEIF(ITY.EQ.13.OR.ITY.EQ.14) THEN                              04220004
         SPHI=1.0                                                       04230004
         EDGE=EINDX2                                                    04240004
         IF(EINDX2.LE.0.02) THEN                                        04250004
           EDGE=FE26ED                                                  04260004
           IF(ITY.EQ.14)EDGE=FE25ED                                     04270004
         ENDIF                                                          04280004
         ENE=ENER                                                       04290004
         IF(ENE.LT.EDGE)RETURN                                          04300004
         IF(ITY.EQ.13)CROSS=6.3E-18/26./26.*(EDGE/ENE)**2.84            04310018
         IF(ITY.EQ.14)CROSS=2.*6.3E-18/25.7/25.7*(EDGE/ENE)**2.84       04320018
C@@      IF(EINDX1.GT.79.)EINDX1=79.                                    04330021
         OPAC=CROSS*10.0**(EINDX1-4.4)                                  04340004
C        OPAC=CROSS*10.0**(EINDX1-4.4)  --- ASSUMES ALLEN               04350004
         SPHI=0.                                                        04360004
         IF(OPAC.LT.100)SPHI=EXP(-OPAC)                                 04370004
         RETURN                                                         04380004
C***********************************  EEKAGEN ********                  04390004
C     ELSEIF(ITY.EQ.13) THEN                                            04400004
C        SPHI=1.                                                        04410004
C        IF(EINDX2.LE.0) EINDX2=FE26ED                                  04420004
C        ENE=ENER*FE26ED/EINDX2                                         04430004
C        IF(ENE.LT.FE26ED)RETURN                                        04440004
C        CROSS=(2./(26.*26.))*(FE26ED/ENE)**2.6666                      04450004
C        OPAC=CROSS*10.0**(EINDX1-17.)                                  04460004
C        SPHI=0.                                                        04470004
C        IF(OPAC.LT.100)SPHI=EXP(-OPAC)                                 04480004
C        RETURN                                                         04490004
C     ELSEIF(ITY.EQ.14) THEN                                            04500004
C        SPHI=1.                                                        04510004
C        IF(EINDX2.LE.0) EINDX2=FE25ED                                  04520004
C        ENE=ENER*FE25ED/EINDX2                                         04530004
C        IF(ENE.LT.FE25ED)RETURN                                        04540004
C        CROSS=(1.3/(26.*26.))*(FE25ED/ENE)**2.6666                     04550004
C        OPAC=CROSS*10.0**(EINDX1-17.)                                  04560004
C        SPHI=0.                                                        04570004
C        IF(OPAC.LT.100)SPHI=EXP(-OPAC)                                 04580004
C        RETURN                                                         04590004
C*************************************************                      04600003
      ELSEIF(ITY.EQ.15) THEN                                            04610000
         SPHI=1.                                                        04620000
         IF(ENER.LE.EINDX1)RETURN                                       04630000
         SPHI=EXP((EINDX1-ENER)/EINDX2)                                 04640021
C@@      IF((EINDX1-ENER)/EINDX2.LT.100)SPHI=EXP((EINDX1-ENER)/EINDX2)  04650021
         RETURN                                                         04660000
      ELSEIF(ITY.EQ.16) THEN                                            04670000
         ENERA=DLOG10(ENER)                                             04680000
         DO 171 I=2,NUMBKG                                              04690000
         IF(ENERA.LT.CBACK(I,2))GO TO 174                               04700000
171      CONTINUE                                                       04710000
         I=NUMBKG                                                       04720000
174      I=I-1                                                          04730000
         SLOPE=(CBACK(I,1)-CBACK(I+1,1))/(CBACK(I,2)-CBACK(I+1,2))      04740000
         TERM2=CBACK(I+1,1)+(ENERA-CBACK(I+1,2))*SLOPE                  04750000
         SPHI=10.**TERM2                                                04760000
CC                                                                      04770000
CC  LINE ABSORPTION -- ITY=20: GAUSSIAN OPACITY (DOPPLER)               04780000
CC                  -- ITY=21: LORENTZIAN OPACITY (DAMPING)             04790000
CC                  -- ITY=22: GAUSSIAN   ABS (TO GET EQUIVALENT WIDTH) 04800000
CC                  -- ITY=23: LORENTZIAN ABS (TO GET EQUIVALENT WIDTH) 04810000
CC                  -- ITY=24: TRIANGULAR ABS (TO GET EQUIVALENT WIDTH) 04820000
CC   DIPOLE= 2.8E-21 KEV*CM**2                                          04830000
CC                                                                      04840000
      ELSEIF(ITY.EQ.20.OR.ITY.EQ.21.) THEN                              04850000
         DE=ENER-EINDX2                                                 04860000
         IF(ITY.EQ.20)THEN                                              04870000
           DE2=DE*DE/EINDX3/EINDX3                                      04880000
           PROFIL=0.                                                    04890000
           IF(DE2.LT.100) PROFIL=EXP(-DE2)/PISQ/EINDX3                  04900000
         ELSEIF(ITY.EQ.21)THEN                                          04910000
           G=EINDX3/4./PI                                               04920000
           PROFIL=G/(DE**2+G**2)/PI                                     04930000
         ENDIF                                                          04940000
         OPAC=DIPOLE*PROFIL*10**(EINDX1-21)                             04950000
         SPHI=0.                                                        04960000
         IF(OPAC.LT.100)SPHI=EXP(-OPAC)                                 04970000
         RETURN                                                         04980000
      ELSEIF(ITY.EQ.22)THEN                                             04990000
        DE2=(ENER-EINDX2)/2./EINDX3                                     05000000
        DE2=DE2*DE2                                                     05010000
        SPHI=0                                                          05020000
        IF(DE2.LT.100) SPHI=EINDX1/1.414213562/PISQ/EINDX3*EXP(-DE2)    05030000
        SPHI=1.-SPHI                                                    05040000
        IF(SPHI.LT.0) SPHI=0.                                           05050000
        RETURN                                                          05060000
      ELSEIF(ITY.EQ.23)THEN                                             05070000
        DE=ENER-EINDX2                                                  05080000
        G=EINDX3/4./PI                                                  05090000
        SPHI=EINDX1/PI*G/(DE*DE+G*G)                                    05100000
        SPHI=1.-SPHI                                                    05110000
        IF(SPHI.LT.0) SPHI=0.                                           05120000
        RETURN                                                          05130000
      ELSEIF(ITY.EQ.24) THEN                                            05140000
        IF(EINDX1.GT.1.) EINDX1=1.                                      05150000
        IF(EINDX1.LT.0.) EINDX1=0.                                      05160000
        DE=ABS(ENER-EINDX2)                                             05170000
        SPHI=1.                                                         05180000
        IF(DE.LT.EINDX3) SPHI=(1.-EINDX1)/EINDX3*DE+EINDX1              05190000
        RETURN                                                          05200000
      ELSEIF(ITY.EQ.30) THEN                                            05210000
C******* FIREMAN OPACITY ***84.1.31   S.K.                              05220001
        CALL FIREMA(ENER,EINDX1,SPHI)                                   05230000
         RETURN                                                         05240000
      ELSEIF(ITY.EQ.31) THEN                                            05250004
C******* USER  OPACITY ***84.2. 8   S.K.                                05260004
        CALL USERAB(ENER,EINDX1,SPHI)                                   05270004
         RETURN                                                         05280004
      ELSEIF(ITY.EQ.32) THEN                                            05290004
C******* A*E**(-ALPH)  ***84.2. 8   S.K.                                05300004
        CALL AEALPH(ENER,EINDX1,EINDX2,EINDX3,SPHI)                     05310004
         RETURN                                                         05320004
      ELSEIF(ITY.EQ.33) THEN                                            05330004
C******* THOMSON SCATTERING **84.2.27   S.K.                            05340004
        CALL THOMS(ENER,EINDX1,SPHI)                                    05350004
         RETURN                                                         05360004
      ELSEIF(ITY.EQ.34) THEN                                            05370004
C******* THOMSON SCATTERING  E DIPENDENT**84.3.12   S.K.                05380004
        CALL THOMS2(ENER,EINDX1,SPHI)                                   05390004
         RETURN                                                         05400004
      ELSEIF(ITY.EQ.35.OR.ITY.EQ.36) THEN                               05410004
C******* COMPTONIZED B-B *********  ******84.8.6  S.K./9.4. S.K.        05420004
C      EINDX1=KT(B.B) KEV (FOR 35)/ ALPH (FOR36 )                       05430004
C      EINDX2=TAU     OPTICAL DEPTH                                     05440000
C      EINDX3=KT(ELECTRON) KEV                                          05450000
C                                                                       05460000
         PROB=1.-EXP(-EINDX2)                                           05470000
         GAMMA=1.0+7.8278E-3*EINDX3+4.5955E-5*EINDX3**2                 05480000
         NPORI=20                                                       05490004
         CALL BCOMP(EINDX1,GAMMA,PROB,ENER,NPORI,FOUT,ITY)              05500004
         SPHI=FOUT                                                      05510000
*      ELSEIF(ITY.EQ.37) THEN                                            05520004
C******* COMPTONIZED P-L *********  ******84.8.6    S.K./9.4 S.K.       05530004
C     COMMON /PARMS/TS,TE,ODEPTH,EF,X                                   05540004
C      EINDX1= ALPH (PHOTON INDEX)                                      05550004
C      EINDX2=TAU     OPTICAL DEPTH                                     05560004
C      EINDX3=KT(ELECTRON) KEV                                          05570004
C                                                                       05580004
*         TS=EINDX1                                                      05590004
*         ODEPTH=EINDX2                                                  05600004
*         TE=EINDX3                                                      05610004
*         EF=ENER                                                        05620004
*         CALL SETK3                                                     05630004
*         CALL CMPSUN(FOUT)                                              05640004
*         SPHI=FOUT                                                      05650004
      ELSEIF(ITY.EQ.38.OR.ITY.EQ.39) THEN                               05660002
C******* COMPTONIZED B-B (OCULTATION)*****84.11.2  S.K./9.4. S.K.       05670002
C      EINDX1=KT(B.B) KEV (FOR 38)/ ALPH (FOR39 )                       05680002
C      EINDX2=TAU     OPTICAL DEPTH                                     05690002
C      EINDX3=KT(ELECTRON) KEV                                          05700002
C                                                                       05710002
         PROB=1.-EXP(-EINDX2)                                           05720002
         GAMMA=1.0+7.8278E-3*EINDX3+4.5955E-5*EINDX3**2                 05730002
         NPORI=20                                                       05740002
         CALL BCOMPO(EINDX1,GAMMA,PROB,ENER,NPORI,FOUT,ITY)             05750002
         SPHI=FOUT                                                      05760002
      ELSEIF(ITY.EQ.40.OR.ITY.EQ.41.OR.ITY.EQ.42) THEN                  05770013
C******* COMPTONIZED B-B (OCULTATION)*****84.11.5    S.K.               05780004
C      EINDX1=KT(B.B) KEV(FOR 40,42)/ ALPH(P.L.) (FOR.41)               05790011
C      EINDX2=TAU     OPTICAL DEPTH                                     05800004
C      EINDX3=KT(ELECTRON) KEV                                          05810004
C       ITY.EQ.42 : OCULTATION                                          05820011
C                                                                       05830004
         PROB=1.-EXP(-EINDX2)                                           05840004
         F1KTMC=1.+2*EINDX3/511.002                                     05850004
         NPORI=20                                                       05860004
         CALL BCOMPR(EINDX1,F1KTMC,PROB,ENER,NPORI,FOUT,ITY)            05870004
         SPHI=FOUT                                                      05880004

      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

      ENDIF                                                             05890004
150   SPHI=SPHI*SCALE                                                   05900000
      RETURN                                                            05910000
      END                                                               05920000
      BLOCK DATA                                                        05930000
      IMPLICIT REAL*8 (A-H,O-Z)                                         05940004
C      COMMON /COEF/ NUMGU,COEFB(7,2),GAUSS(5,2),BREAK(7),NUMBKG,        05950000
C     *               CBACK(15,2)                                        05960000
      COMMON /COEF/ COEFB(7,2),GAUSS(5,2),BREAK(7), CBACK(15,2)
      COMMON /NCOEF/ NUMGU,NUMBKG

      DATA NUMGU /5/                                                    05970000
      DATA COEFB /2.,1.9066,1.74,1.59,1.47,1.30,1.,                     05980000
     *            .7,.74666,0.88,1.08,1.32,2.00,5./                     05990000
      DATA GAUSS /.236926885056189,.478628670499366,.568888888888888,   06000000
     *            .478628670499366,.236926885056189,                    06010000
     *           -.906179845938664,-.538469310105683,0.0,               06020000
     *            .538469310105683, .906179845938664/                   06030000
      DATA BREAK /1.E38,2.,1.25,.75,.5,.25,.1/                          06040000
      DATA NUMBKG/7/                                                    06050000
      DATA CBACK /.5,-.5,-1.5,-2.5,-3.85,-4.9,-5.65,8*0.,               06060000
     *            0.,1.,1.6,2.,2.47,3.,3.47,8*0./                       06070000
      END                                                               06080000
      FUNCTION SECTON(ENER,COEFAT)                                      06090000
      IMPLICIT REAL*8 (A-H,O-Z)                                         06100004
C                                                                       06110000
C     THE PURPOSE OF SECTON IS TO EVALUATE AN CROSS SECTION FORMULA     06120000
C                                                                       06130000
      DIMENSION COEFAT(4)                                               06140000
      ELNB=DLOG(ENER)                                                   06150000
      PSUM=0.                                                           06160000
      DO 118 I=1,4                                                      06170000
      ELN=ELNB**(I-1)                                                   06180000
      PSUM=PSUM+ELN*COEFAT(I)                                           06190000
118   CONTINUE                                                          06200000
      SECTON=EXP(PSUM)                                                  06210000
      RETURN                                                            06220000
      END                                                               06230000
C********************************************************************** 06240000
C        CHAPLINE & STEVENS (A.P.J. 184 1041(1973))                     06250000
C     COMPUTON SCATTERING NI YOTE BREMSSTRAHRUNG KARA ZURERU            06260000
C        KINGI SHIKI                                                    06270000
C          ELECTRON TEMPERATURE  ;T                                     06280000
C          THOMSON SCATTERING LENGTH  ;R                                06290000
C          PHOTON ENERGY  ;  E                                          06300000
C********************************************************************** 06310000
      SUBROUTINE CHASTE(E,T,R,CR)                                       06320000
      IMPLICIT REAL*8 (A-H,O-Z)                                         06330004
C---------------------------------------------------------------------- 06340000
C     INPUT ARGUMENTS                                                   06350000
C        E  :  PHOTON ENERGY (KEV)                                      06360000
C        T  :  ELECTRON TENPERATURE (KEV)                               06370000
C        R  :  THOMSON SCATTERING LENGTH                                06380000
C     OUTPUT AREGUMENT                                                  06390000
C        CR :  CRECTION FACTER OF FREE FREE RADIATION                   06400000
C                                                                       06410000
C        CR=1+B-C+D                                                     06420000
C            A=2*T/510                                                  06430000
C            ARS=A*R**2                                                 06440000
C            EKT=E/T ,EKT2=(E/T)**2 ,EKT3=(E/T)**3                      06450000
C            P=3*(EKT3-EKT2-EKT)/(2*ARS)                                06460000
C            Q=3*(2*EKT2-EKT3)/(4*ARS)                                  06470000
C           B=3*EKT3*ARS/8                                              06480000
C           C=P*(EXP(-ARS)-1+ARS)                                       06490000
C           D=Q*(EXP(-ARS)*(ARS+1)-1)                                   06500000
C-----------------------------------------------------------------      06510000
             A=2.0*T/510.0                                              06520000
             ARS=A*R**2.0                                               06530000
              IF(ARS.LT.0.000002) THEN                                  06540000
                CR=1.0                                                  06550000
                RETURN                                                  06560000
              ENDIF                                                     06570000
             EKT=E/T                                                    06580000
             EKT2=(E/T)**2                                              06590000
             EKT3=(E/T)**3                                              06600000
             P=3.0*(EKT3-EKT2-EKT)/(2.0*ARS)                            06610000
             Q=3.0*(2.0*EKT2-EKT3)/(4.0*ARS)                            06620000
            B=3.0*EKT3*ARS/8.0                                          06630000
            C=P*(EXP(-ARS)-1.0+ARS)                                     06640000
            D=Q*(EXP(-ARS)*(ARS+1.0)-1.0)                               06650000
         CR=1+B-C+D                                                     06660000
         RETURN                                                         06670000
         END                                                            06680000
C**********************************************************             06690004
        SUBROUTINE FIREMA(ENER,EINDX1,SPHI)                             06700000
      IMPLICIT REAL*8 (A-H,O-Z)                                         06710004
C---------------------------------------------------------------------  06720000
C        FIREMAN NO CROSS SECTION O TSUKAU  (FIREMAN 1974,AP.J.,187,57  06730000
C           GAS NO TOKI                                                 06740000
C            ELEMENT  O,NE,MG,SI,S,AR,CA,CR,FE,NI                       06750003
C                                                                       06760003
C------------------------------------------------------------------     06770000
         DIMENSION EDGE(11),COEF(11),EIND(11)                           06780002
C----------- K-EDGE ---------------------------------------------       06790002
         DATA EDGE/  0.30,  0.53,  0.87,  1.30,  1.84,  2.47,           06800002
     &               3.20,  4.04,  5.99,  7.11,  8.33/                  06810002
C----------- INDEX ----------------------------------------------       06820002
         DATA EIND/ -2.82, -2.90, -2.90, -2.85, -2.70, -2.90,           06830004
     &              -2.80, -2.81, -2.72, -2.65, -2.76/                  06840004
C----------- KEISU *10**23 --------------------------------------       06850002
         DATA COEF/  9.53, 20.25, 24.34, 26.44, 26.65, 35.62,           06860002
     &              34.08, 35.28, 30.59, 52.85, 70.34/                  06870002
C-----------------------------------------------------------------      06880000
          IF(ENER.LT.EDGE(1)) THEN                                      06890000
           CROSS=2.11*ENER**(-2.6)                                      06900000
           GO TO 3000                                                   06910000
          ENDIF                                                         06920000
 1001     DO 1000 I1=1,10                                               06930003
          I1A=I1                                                        06940000
         IF(ENER.GE.EDGE(I1).AND.ENER.LT.EDGE(I1+1))GO TO 2000          06950000
 1000     CONTINUE                                                      06960000
          I1A=11                                                        06970003
 2000     CONTINUE                                                      06980000
            CROSS=COEF(I1A)*ENER**(EIND(I1A))                           06990002
 3000     OPAC=10.0**(EINDX1-23)*CROSS                                  07000000
           SPHI=0.0                                                     07010000
         IF(OPAC.LT.100)SPHI=EXP(-OPAC)                                 07020000
         RETURN                                                         07030000
         END                                                            07040000
C**********************************************************             07050004
        SUBROUTINE USERAB(ENER,EINDX1,SPHI)                             07060004
      IMPLICIT REAL*8 (A-H,O-Z)                                         07070004
C---------------------------------------------------------------------  07080004
C        USER  NO CROSS SECTION DE KYUSHU KEISU O TSUKURU               07090004
C                                                                       07100004
C          ** 84/12/14 .... FOR  PARTIAL IONIZED GASS                   07110023
C           MG,SI,S,AR,CA,CR,FE (HE LIKE)                               07120025
C              FIREMAN TO INDEX=2.84                                    07121025
C------------------------------------------------------------------     07130004
         DIMENSION EDGE(11),COEF(11),EIND(11)                           07140004
C----------- K-EDGE ---------------------------------------------       07150023
CHE-LIKE DATA EDGE/  1.76,  2.44,  3.08,                                07160027
C    &               4.12,  5.13,  7.48,  8.83, 99.00,                  07170027
C    &               99.0,  99.0,  99.0/                                07180027
C----------- K-EDGE ---------------------------------------------       07181027
         DATA EDGE/  1.96,  2.67,  3.50,                                07182027
     &               4.42,  5.47,  7.90,  8.83, 99.00,                  07183027
     &               99.0,  99.0,  99.0/                                07184027
C----------- INDEX ----------------------------------------------       07190023
         DATA EIND/ -2.84, -2.84, -2.84,                                07200025
     &              -2.84, -2.84, -2.84, -2.84, -2.84,                  07210025
     &              -2.84, -2.84,  99.0/                                07220025
C----------- KEISU *10**23 --------------------------------------       07230023
         DATA COEF/  2.28,  6.60, 11.61,                                07261025
     &              14.40, 16.20, 18.50, 56.50, 99.00,                  07262025
     &               99.0,  99.0,  99.0/                                07263024
C-----------------------------------------------------------------      07270023
          IF(ENER.LT.EDGE(1)) THEN                                      07280004
           CROSS=0.00                                                   07290023
           GO TO 3000                                                   07300004
          ENDIF                                                         07310004
 1001     DO 1000 I1=1,10                                               07320004
          I1A=I1                                                        07330004
         IF(ENER.GE.EDGE(I1).AND.ENER.LT.EDGE(I1+1))GO TO 2000          07340004
 1000     CONTINUE                                                      07350004
          I1A=11                                                        07360004
 2000     CONTINUE                                                      07370004
            CROSS=COEF(I1A)*ENER**(EIND(I1A))                           07380004
 3000     OPAC=10.0**(EINDX1-23)*CROSS                                  07390004
           SPHI=0.0                                                     07400004
         IF(OPAC.LT.100)SPHI=EXP(-OPAC)                                 07410004
         RETURN                                                         07420004
         END                                                            07430004
C****************************************************                   07440004
       SUBROUTINE AEALPH(ENER,EINDX1,EINDX2,EINDX3,SPHI)                07450004
      IMPLICIT REAL*8 (A-H,O-Z)                                         07460004
C---------------------------------------------------                    07470004
C        OPECITY O A*E**(-ALPH) TO SURU                                 07480004
C               EINDX1...... NH                                         07490004
C               EINDX2...... A                                          07500004
C               EINDX3...... ALPH                                       07510004
C----------------------------------------------------                   07520004
          CROSS=EINDX2*ENER**(-EINDX3)                                  07530004
         OPAC=10.0**(EINDX1-22)*CROSS                                   07540004
         SPHI=0.                                                        07550004
         IF(OPAC.LT.100)SPHI=EXP(-OPAC)                                 07560004
         RETURN                                                         07570004
         END                                                            07580004
C*****************************************************                  07590004
        SUBROUTINE THOMS(ENER,EINDX1,SPHI)                              07600004
      IMPLICIT REAL*8 (A-H,O-Z)                                         07610004
C-----------------------------------------------------                  07620004
C        NH ATARI NH NO ELECTRON                                        07630004
C             ELECTRON 1 KO ATARI .... 6.65*10**(-25) CM**2             07640004
C             ELECTRON  10**22 ATARI...6.65*10**(-3)  CM**2             07650004
C                                                                       07660004
C------------------------------------------------------                 07670004
        DATA THO22 /0.00665/                                            07680004
         OPAC=10**(EINDX1-22)*THO22                                     07690004
C        OPAC=10**(EINDX1-22)*THO22*1.2                                 07700004
         SPHI=0.                                                        07710004
         IF(OPAC.LT.100)SPHI=EXP(-OPAC)                                 07720004
         RETURN                                                         07730004
         END                                                            07740004
C*****************************************************                  07750004
        SUBROUTINE THOMS2(ENER,EINDX1,SPHI)                             07760004
      IMPLICIT REAL*8 (A-H,O-Z)                                         07770004
C-----------------------------------------------------                  07780004
C        NH ATARI NH NO ELECTRON                                        07790004
C         PHAI0=  ELECTRON 1 KO ATARI .... 6.65*10**(-25) CM**2         07800004
C              =  ELECTRON  10**22 ATARI...6.65*10**(-3)  CM**2         07810004
C        W.HEITLER P.226                                                07820004
C            PHAI=PHAI0*(1-2.0*GAMMA+5.2*GAMMA**2)                      07830004
C                 GAMMA=ENER(KEV)/512(KEV)                              07840004
C------------------------------------------------------                 07850004
        DATA THO22 /0.00665/                                            07860004
        DATA ELEC0 /512.0/                                              07870004
         GAMMA=ENER/ELEC0                                               07880004
         OPAC=10**(EINDX1-22)*THO22*(1.0-2.0*GAMMA+5.2*GAMMA**2)        07890004
         SPHI=0.                                                        07900004
         IF(OPAC.LT.100)SPHI=EXP(-OPAC)                                 07910004
         RETURN                                                         07920004
         END                                                            07930004
C***************************************************************        07940000
       SUBROUTINE BCOMP(T,G,P,E,N,FOUT,ITY)                             07950004
       IMPLICIT REAL*8 (A-H,O-Z)                                        07960002
C--------------------------------------------------------------         07970000
C      COMPTONIZED B-B NO SPECTRUM O MOTOMERU                           07980000
C       INPUT AREGUMENTS                                                07990000
C          T  : B-B TEMP.(KEV)FOR 35/ ALPH FOR 36                       08000004
C          G  : (1+4*KTE/M/C**2+12*(KTE/M/C**2)**2) KTE: ELECTRON TE-   08010000
C               1.0+7.8278E-3*KTE+4.5955E-5*KTE**2    MPERATURE(KEV)    08020000
C          P  : PROBABIRITY OF SCATTERING (1-EXP(-TAU) TAU:OPTICAL      08030000
C                                                          DEPTH)       08040000
C          E  : HOSHI ENERGY (KEV)                                      08050000
C          N  : PORINOMIAL NO KAISUU                                    08060000
C          ITY: INPUT SPECTRUM TYPE  35 : BB/ 36 : PL                   08070004
C       OUTPUT AREGUMENTS                                               08080000
C          FOUT : ENERGY E NO PHOTON SUU                                08090000
C                                                                       08100000
C        TAU<1.0 DE NAITO TSUKAENAI                                     08110000
C                                        N.SHIBAZAKI & K.MITSUDA        08120000
C                              SUBROUTINE: KITAMOTO                     08130000
C-------------------------------------------------------------          08140000
        BB(T,EGN)=EGN*EGN/(EXP(EGN/T)-1.)                               08150000
        PL(T,EGN)=EGN**(-T)                                             08160004
C--------------- TAU CHECK -----------------                            08170000
         IF(P.GT.1.0)GO TO 999                                          08180000
C--------------- ITY CHECK -----------------                            08190004
         IF(ITY.NE.35.AND.ITY.NE.36)GO TO 998                           08200004
C--------------------------------------------                           08210000
   1    CONTINUE                                                        08220004
        FOUT=0.                                                         08230000
 1001   DO 1000 I1=0,N                                                  08240000
         FI1=I1                                                         08250001
         EGN=E/G**FI1                                                   08260001
         IF(EGN.LT.0.000001)GO TO 9000                                  08270004
         IF(ITY.EQ.35)FIN=BB(T,EGN)                                     08280004
         IF(ITY.EQ.36)FIN=PL(T,EGN)                                     08290004
         FOUT=FOUT+FIN*(P/G)**FI1*(1.-P)                                08300001
 1000   CONTINUE                                                        08310000
 9000   RETURN                                                          08320004
C--------------- MADA DEKITENAI -----------------                       08330004
 998      PRINT *,'??????? ITY = ',ITY,'  IN BCOMP ?????'               08340004
           STOP                                                         08350004
C--------------- TAU G.T. 1 ---------------------                       08360004
 999      PRINT *,'??????? TAU G.T.1 IN BCOMP ??????????'               08370004
          PRINT *,'??????? SET TAU.EQ.1 BCOMP ??????????'               08380000
C         P=1.0                                                         08390004
          GO TO 1                                                       08400004
C------------------------------------------------                       08410000
        END                                                             08420000
C***************************************************************        08430002
       SUBROUTINE BCOMPO(T,G,P,E,N,FOUT,ITY)                            08440002
       IMPLICIT REAL*8 (A-H,O-Z)                                        08450002
C--------------------------------------------------------------         08460002
C      COMPTONIZED B-B NO SPECTRUM O MOTOMERU                           08470002
C       INPUT AREGUMENTS                                                08480002
C          T  : B-B TEMP.(KEV)FOR 35/ ALPH FOR 36                       08490002
C          G  : (1+4*KTE/M/C**2+12*(KTE/M/C**2)**2) KTE: ELECTRON TE-   08500002
C               1.0+7.8278E-3*KTE+4.5955E-5*KTE**2    MPERATURE(KEV)    08510002
C          P  : PROBABIRITY OF SCATTERING (1-EXP(-TAU) TAU:OPTICAL      08520002
C                                                          DEPTH)       08530002
C          E  : HOSHI ENERGY (KEV)                                      08540002
C          N  : PORINOMIAL NO KAISUU                                    08550002
C          ITY: INPUT SPECTRUM TYPE  35 : BB/ 36 : PL                   08560002
C       OUTPUT AREGUMENTS                                               08570002
C          FOUT : ENERGY E NO PHOTON SUU                                08580002
C                                                                       08590002
C        TAU<1.0 DE NAITO TSUKAENAI                                     08600002
C                                        N.SHIBAZAKI & K.MITSUDA        08610002
C                              SUBROUTINE: KITAMOTO                     08620002
C-------------------------------------------------------------          08630002
        BB(T,EGN)=EGN*EGN/(EXP(EGN/T)-1.)                               08640002
        PL(T,EGN)=EGN**(-T)                                             08650002
C--------------- TAU CHECK -----------------                            08660002
         IF(P.GT.1.0)GO TO 999                                          08670002
C--------------- ITY CHECK -----------------                            08680002
         IF(ITY.NE.38.AND.ITY.NE.39)GO TO 998                           08690002
C--------------------------------------------                           08700002
   1    CONTINUE                                                        08710002
        FOUT=0.                                                         08720002
 1001   DO 1000 I1=1,N                                                  08730002
         FI1=I1                                                         08740002
         EGN=E/G**FI1                                                   08750002
         IF(EGN.LT.0.000001)GO TO 9000                                  08760002
         IF(ITY.EQ.38)FIN=BB(T,EGN)                                     08770002
         IF(ITY.EQ.39)FIN=PL(T,EGN)                                     08780002
         FOUT=FOUT+FIN*(P/G)**FI1*(1.-P)                                08790002
 1000   CONTINUE                                                        08800002
 9000   RETURN                                                          08810002
C--------------- MADA DEKITENAI -----------------                       08820002
 998      PRINT *,'??????? ITY = ',ITY,'  IN BCOMPO?????'               08830002
           STOP                                                         08840002
C--------------- TAU G.T. 1 ---------------------                       08850002
 999      PRINT *,'??????? TAU G.T.1 IN BCOMPO??????????'               08860002
          PRINT *,'??????? SET TAU.EQ.1 BCOMPO??????????'               08870002
C         P=1.0                                                         08880002
          GO TO 1                                                       08890002
C------------------------------------------------                       08900002
        END                                                             08910002
C***************************************************************        08920004
       SUBROUTINE BCOMPR(T,G,P,E,N,FOUT,ITY)                            08930004
       IMPLICIT REAL*8 (A-H,O-Z)                                        08940004
C--------------------------------------------------------------         08950004
C      COMPTONIZED B-B NO SPECTRUM O MOTOMERU                           08960004
C       INPUT AREGUMENTS                                                08970004
C          T  : B-B TEMP.(KEV)FOR 35/ ALPH FOR 36                       08980004
C          G  : 1.+2*KTE/M/C**2    KTE: ELECTRON TE                     08990004
C          P  : PROBABIRITY OF SCATTERING (1-EXP(-TAU) TAU:OPTICAL      09000004
C                                                          DEPTH)       09010004
C          E  : HOSHI ENERGY (KEV)                                      09020004
C          N  : PORINOMIAL NO KAISUU                                    09030004
C          ITY: INPUT SPECTRUM TYPE  35 : BB/ 36 : PL                   09040004
C       OUTPUT AREGUMENTS                                               09050004
C          FOUT : ENERGY E NO PHOTON SUU                                09060004
C                                                                       09070004
C        TAU<1.0 DE NAITO TSUKAENAI                                     09080004
C                    MODIFIED N.SHIBAZAKI & K.MITSUDA                   09090004
C                              SUBROUTINE: KITAMOTO                     09100004
C-------------------------------------------------------------          09110004
        DATA MC2/511.002/                                               09120004
        BB(T,EGN)=EGN*EGN/(EXP(EGN/T)-1.)                               09130004
        PL(T,EGN)=EGN**(-T)                                             09140004
C--------------- TAU CHECK -----------------                            09150004
         IF(P.GT.1.0)GO TO 999                                          09160004
C--------------- ITY CHECK -----------------                            09170004
         IF(ITY.NE.40.AND.ITY.NE.41.AND.ITY.NE.42)GO TO 998             09180012
C--------------------------------------------                           09190004
C         PRINT *,' INPUT PARAMETERS IN BCOMPR '                        09200010
C         PRINT *,'  T,G,P,E,N,ITY ',T,G,P,E,N,ITY                      09210010
   1    CONTINUE                                                        09220004
        IF(ITY.EQ.40)FIN=BB(T,E)                                        09230004
        IF(ITY.EQ.41)FIN=PL(T,E)                                        09240004
        EB=E                                                            09250004
        DE=1.                                                           09260004
        FOUT=FIN*(1.-P)                                                 09270004
        IF(ITY.EQ.42)FOUT=0                                             09280011
C        PRINT *,' FOUT 0 ',FOUT                                        09290010
 1001   DO 1000 I1=1,N                                                  09300004
         FI1=I1                                                         09310004
         GEBMC2=G-EB/MC2                                                09320004
         EBB=EB/G/GEBMC2                                                09330004
         DE=DE/GEBMC2**2                                                09340005
         PFI=P**FI1                                                     09350007
C          PRINT *,'  I1,GEBMC2,EBB,DE,PFI ',I1,GEBMC2,EBB,DE,PFI       09360010
         IF(DE.LT.0.000001)GO TO 9000                                   09370004
         IF(PFI.LT.0.000001)GO TO 9000                                  09380007
         IF(EBB/T.GT.168.)EBB=168.*T                                    09390017
         IF(ITY.EQ.40.OR.ITY.EQ.42)FIN=BB(T,EBB)                        09400014
         IF(ITY.EQ.41)FIN=PL(T,EBB)                                     09410004
         FOUT=FOUT+FIN*PFI*(1.-P)*DE                                    09420007
C          PRINT *,' FOUT ',FOUT                                        09430010
         EB=EBB                                                         09440004
 1000   CONTINUE                                                        09450004
 9000   RETURN                                                          09460004
C--------------- MADA DEKITENAI -----------------                       09470004
 998      PRINT *,'??????? ITY = ',ITY,'  IN BCOMPR ?????'              09480004
           STOP                                                         09490004
C--------------- TAU G.T. 1 ---------------------                       09500004
 999      PRINT *,'??????? TAU G.T.1 IN BCOMPR ??????????'              09510004
          GO TO 1                                                       09520004
C------------------------------------------------                       09530004
        END                                                             09540004
C                                                                       09550026
C                                                                       09560026
C                                                                       09570026
C**********************************************************             
C  DISK MODEL BY K.MITSUDA(PBSG051.BBSMODEL.FORT77(DISKSPEC)            
C                                                                       
      SUBROUTINE DSKSPC(T,E,FLUX)                                       
      REAL*8    T,E,FLUX                                                
                                                                        
C   SPECTRUM OF ACCRETION DISK                                          
C     T:  TEMPERATURE OF INNER EDGE                                     
C     E:  PHOTON ENERGY                                                 
C     F:  FLUX    (MUST BE NORMALIZED BY SOME CONSTANT)                 
C  norm= (8 pi rin^2 cso i) /(3 d^2) 
                                                                        
      COMMON /BBFLD/ ET                                                 
      REAL*8    ET                                                      
      EXTERNAL  BBFOLD                                                  
      REAL*8    BBFOLD,A/0.0/,B/1.0/,EPSA/0.0/,EPSR/0.003/,ERR          
      INTEGER*4 NMIN/20/,NMAX/641/,N,ICON                               
                                                                        
      ET=E/T
*      print *,"ET=",ET
*      CALL DAQE(A,B,BBFOLD,EPSA,EPSR,NMIN,NMAX,FLUX,ERR,N,ICON)        
      call bbint(a,b,flux,icon)

*      print *,"a,b,flux=",a,b,flux

      IF(ICON.NE.0) THEN                                                
        WRITE(6,*) 'bbint ERROR: ',ICON                                  
      ENDIF                                                             
      IF(ICON.GE.30000) THEN                                            
        WRITE(6,*) 'bbint ERROR: ',ICON                                  
        STOP                                                            
      ENDIF                                                             
      FLUX=FLUX*E*E                                                     
                                                                        
      RETURN                                                            
      END                                                               


      subroutine bbint(a,b,flux,icon)
      IMPLICIT REAL*8 (A-H,O-Z)                                         
      integer*4 icon

      COMMON /COEF/ COEFB(7,2),GAUSS(5,2),BREAK(7), CBACK(15,2)
      COMMON /NCOEF/ NUMGU,NUMBKG

      nint=10
      DB=(b-a)/nint
      DDB=DB/2.

*      print *,"db=",db," ddb=",ddb
      SUM=0.                
      DO 100 I=1,nint
         APBB=(2*I-1)*DDB+a
         DO 100 II=1,NUMGU         
         E=DDB*GAUSS(II,2)+APBB
         term=bbfold(E)
         SUM=SUM+GAUSS(II,1)*TERM 
*         print *,"sum=",sum," nint=",nint," ii=",ii

 100  CONTINUE           
*      flux = sum  * bug found on 2004.9.8
      flux = sum*DDB

      icon=0
      return
      end

      FUNCTION BBFOLD(X)                                                
      REAL*8   BBFOLD,X                                                 
                                                                        
      COMMON /BBFLD/ ET                                                 
      REAL*8   ET                                                       
                                                                        
      REAL*8   BEKI,W                                                   
      PARAMETER (BEKI=-11.0D0/3.0D0)                                    
                                                                        
      W=ET/X                                                            
      IF(W.GT.170.0) THEN                                               
        BBFOLD=X**BEKI*EXP(-W)                                          
       ELSE                                                             
        BBFOLD=X**BEKI/(EXP(W)-1.0D0)                                   
      ENDIF                                                             
                                                                        
      RETURN                                                            
      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
