c$$$      implicit none
c$$$      INTEGER NE
c$$$      parameter (ne = 1000)
c$$$      REAL    EAR(0:NE), PARAM(5), PHOTAR(NE), worksp(1)
c$$$      integer i
c$$$
c$$$      do i = 1, ne
c$$$         ear(i) = real(0.01*i)
c$$$         photar(i) = 0.0
c$$$      end do
c$$$
c$$$      param(1) = 10.0
c$$$      param(2) = 2.0
c$$$      param(3) = 0.1
c$$$      param(4) = 50.0
c$$$      param(5) = 50.0
c$$$
c$$$      call pcygni(EAR,NE,PARAM,WORKSP,PHOTAR)
c$$$
c$$$      do i = 1, ne
c$$$         write(*,*) ear(i), photar(i)
c$$$      end do
c$$$
c$$$      end 
c$$$
      SUBROUTINE pcygni(EAR,NE,PARAM,WORKSP,PHOTAR)
      implicit none
      INTEGER NE
      REAL    EAR(0:NE), PARAM(5), WORKSP(*), PHOTAR(NE)

c     param(1) --- continuum level (in photons/keV)
c     param(2) --- E0 (keV)
c     param(3) --- velocity of the flow (in v/c)
c     param(4) --- emission equivalent width (in eV)
c     param(5) --- absorption equivalent width (in eV)

      integer i

      real cont, E0, vel, Eeq, Aeq, Ecenter

      cont  = param(1)
      E0    = param(2)
      vel   = param(3)
      Eeq   = param(4)/1000.0
      Aeq   = param(5)/1000.0

      do i = 1, ne
         Ecenter = (ear(i-1)+ear(i)) / 2.0
         if(Ecenter.gt.E0*(1.0-vel).and.
     $      Ecenter.lt.E0*(1.0+vel)) then
            photar(i) = (cont + Eeq*cont/(2.0*vel*E0))
     $           *(ear(i)-ear(i-1))
         elseif(Ecenter.gt.E0*(1.0+vel).and.
     $          Ecenter.lt.E0*(1.0+vel)+Aeq) then
            photar(i) = 0.0 
         else
            photar(i) = param(1)*(ear(i)-ear(i-1))
         endif
      end do
      END