      subroutine xsismheh(ear,ne,param,Ifl,photar) 
c
c
c     number of model parameters: 8  
c     1: Te, electron temperature in keV  
c        if Te > 0 exact redistribution function (RF),
c        if Te < 0 isotropic approximation for RF
c     2: Tbb, blackbody temperature at in the inner edge in keV  
c     3: tau, vertical optical depth of the corona
c        if tau > 0 isotropic source function (SF) starting
c                   from the 3rd scattering
c        if tau < 0 isotropic SF starting from the 1st scattering
c     4: cosIncl, cosine of inclination angle
c     5: Rin, inner radius of the disk 
c     6, R, amount of reflection (if < 0 then only reflection component) 
c     7: FeAb, iron abundance in units of solar
c     8: redshift 
c
c     algorithm:
c           exact solution of the radiative transfer equation
c           in HEMISPHERE geometry, 
c                     neglecting the energy balance 
c electrons: Maxwellian 
c           includes: Comptonization of the BB radiation from the 
c           disk, and reflection
c           
      implicit none
      integer np,i,nem
      PARAMETER(nem=5000) 
      PARAMETER(np=20)
      integer ne,iigeom,ne1,iidistf,iion,Ifl 
      real*4  param(*),ear(0:ne),photar(ne)
      real*4  pa0(np)
      real*8  pac(np)    
      real*8  phcon(nem),phref(nem),phblb(nem),obj(nem)  
      real*4  phtot(nem)
      logical firstcall
      save firstcall
      data firstcall/.true./
c      logical fl0
      data    pa0/np*9999./

      if (firstcall) then
        call xwrite('ISM Version 2.3',5) 
        call xwrite('Comptonization in hemisphere-corona',5)
        call xwrite('sombrero model',5)
        call xwrite('Maxwellian electrons',5)
        call xwrite('Questions: Juri Poutanen',5)
        call xwrite('juri@astro.su.se',5)
        call xwrite('Attention ! This code is under development',5)
        call xwrite('and is NOT sufficiently tested yet.',5)
        call xwrite('Before publishing any results ask the author',5)
        call xwrite('for the latest version of the code.',5)
        call xwrite('Your comments are welcome.',5)
        firstcall = .false.
      end if

      ne1=ne+1 
c check if the parameters were changed  
c      do n=1,np
c         if(param(n).ne.pa0(n)) fl0=.true.
c      enddo
c electron temperature
         pac(1)=dble(param(1))
c power-law index
         pac(2)=0d0
c gamma_min
         pac(3)=0d0
c gamma_max
         pac(4)=0d0
c disc temperature (keV)
         pac(5)=dble(param(2))
c Thomson optical depth
         pac(6)=dble(param(3))
c cosine of inclination angle
         pac(7)=dble(param(4))
c g, or Rin (inner radius of the disc)
         pac(8)=dble(param(5))
c normalization of reflection
         pac(9)=dabs(dble(param(6)))
c iron abundance
         pac(10)=dble(param(7))
c "metal" abundances
         pac(11)=0d0
c xi (ionization parameter)
         pac(12)=0d0
c disc temp (for ionization)
         pac(13)=0d0
c redshift
         pac(14)=dble(param(8))
         pac(15)=0d0
         pac(16)=0d0
         pac(17)=0d0
         pac(18)=0d0
         pac(19)=0d0
         pac(20)=0d0



C IGEOM=0, 3  SLAB geometry
C IGEOM=1, 4 CYLINDER geometry (with T0TAUG=H/R=1/2)
C IGEOM=2, 5 HEMISPHERE geometry  (TOTAUG=1, H=R)
c for IIGEOM.gt.2 -- holly disk
c HEMISPHERE
      IIGEOM=5
C IDISTF=0 Maxwellian electron distribution
C IDISTF=1 monoenergetic electron distribution
C IDISTF=2 pow-law  electron distribution
      IIDISTF=0
c IION=0 neutral reflector (H and He are always ionized)
c IION=1 ionized reflector
      IION=0

c object's frame energies (in mc2) 
         do i=0,ne
            obj(i+1)=dble((1.+param(8))*ear(i)/511.) 
         enddo 
c ---- 
      call ismco(pac,obj,phcon,phref,phblb,ne1,iigeom,iidistf,iion)
c ----
c keep the new parameters  
c         do n=1,np
c            pa0(n)=param(n)
c         enddo
c         fl0=.false.
c      endif
      do i=1,ne1
      if(param(6).ge.0.) then 
c       phtot(i)=sngl(phcon(i)+phref(i)+phblb(i)) 
       phtot(i)=sngl(phcon(i)+phref(i)) 
      else 
       phtot(i)=sngl(phref(i)) 
      endif 
      enddo 
c integrate over the bin size 
      do i=1,ne
        photar(i)=(phtot(i+1)+phtot(i))*(ear(i)-ear(i-1))/2. 
      enddo
 

      return
      end
 
c ------------------------------------------------------------------- c
