
      x = abs( dv ) / dopwid
      y = lorwid / dopwid

      t = cmplx( y, -x )
      s = abs( x ) + y


      if( s .ge. 15.0d0 ) then
         !
         ! region I
         !
         w4 = t * 0.5641896d0 / ( 0.5d0 + t * t )
      else if( s .ge. 5.5d0 ) then
         !
         ! region II
         !
         u  = t * t
         w4 = t * ( 1.410474d0 + u * 0.5641896d0 ) &
              / ( 0.75d0 + u * ( 3.0d0 + u ) )
      else if( y .ge. 0.195d0 * abs( x ) - 0.176d0 ) then
         !
         ! region III
         !
         w4 = ( 16.4955d0 &
              + t * ( 20.20933d0 + t * ( 11.96482d0 &
              + t * ( 3.778987d0 + t * 0.5642236d0 ) ) ) ) &
              / ( 16.4955d0 &
              + t * ( 38.82363d0 + t * ( 39.27121d0 &
              + t * ( 21.69274d0 + t * ( 6.699398d0 + t ) ) ) ) )
      else
         !
         ! region IV
         !
         u  = t * t
!!$         w4 = cexp( u ) &
         w4 = exp( u ) &
              - t * ( 36183.31d0 &
              - u * ( 3321.9905d0 - u * ( 1540.787d0 &
              - u * ( 219.0313d0  - u * ( 35.76683d0 &
              - u * ( 1.320522d0 - u * 0.56419d0 ) ) ) ) ) ) &
              / ( 32066.6d0 &
              - u * ( 24322.84d0 - u * ( 9022.228d0 &
              - u * ( 2186.181d0 - u * ( 364.2191d0 &
              - u * ( 61.57037d0 - u * ( 1.841439d0 - u ) ) ) ) ) ) )
      end if


      lsf = real( w4 )
!!$      lsf = real( w4 ) / ( dopwid * sqrt( pi ) )
