!c===============================================================
      subroutine ibis_eresp(filename,crval1,crval2,crpix1,crpix2,cdelt1,cdelt2,crota2,erange,image,skyexp0,nxm,nx,ny,fast,status)
!c--------------------------------------------------------------         
!c     write image to FITS file
!c
!c     Created: Mon May 16 14:21:02 EDT 1994
!c-------------------------------------------------------------- 
      implicit none
      character*(*) filename
      integer nxm,nx,ny,erange,n,m
       logical anyf

      integer lnblnk,fast
      integer bitpix,pcount
      real*8 norma0,norma1,slope0,slope1,dist,e1,e2,texp
      real*8 crval1,crval2,crpix1,crpix2,cdelt1,cdelt2,crota2
      integer naxes(3),naxis,gcount,extend
      integer unit,status,block,extver,hdutype,simple
      character*400 comment
      integer nmapm,nmapx,nmapy
      parameter(nmapm=2000)

      real skyexp(nmapm,nmapm)
      real skyexp0(nmapm,nmapm)
      real image(nmapm,nmapm)

      status=0
      norma1=0.0
      slope1=1.0
      if(erange .eq. -1) then
         norma0=0.0
         slope0=1.0
      else if(erange .eq. 1) then
         norma0=0.48859433
         slope0=1.3304603
         norma1=0.083260782
         slope1=0.68063210
      else if(erange .eq. 2) then
         norma0=0.21444822
         slope0=1.7809195
         norma1=0.045216999
         slope1=1.3624652
      else if(erange .eq. 3) then
         norma0=0.097560597
         slope0=1.0071800
      else if(erange .eq. 9) then
         norma0=0.32660987
         slope0=1.7066450
         norma1=0.24734285
         slope1=0.17543571
      else if(erange .eq. 10) then
         norma0=0.10310040
         slope0=1.3237412
      else if(erange .eq. 8) then
         norma0=-0.079264878
         slope0=0.91217010
      else if(erange .eq. 13) then
         norma0=-0.035629722
         slope0=1.3605498
      else if(erange .ge. 4 .and. erange .le. 7) then
         norma0=0.0
         slope0=1.0
      else if(erange .eq. 11 .or. erange .eq. 12 .or. erange .eq. 14 .or. erange .eq. 15) then
         norma0=0.0
         slope0=1.0
      else
         print*,'No such energy interval ',erange
         stop
      endif

      call ftgiou(unit,status)
      call ftopen(unit,filename,0,block,status)
      if(status.ne.0) then
         print*,'Bad file, skip it'
         call ftclos(unit,status)
         call ftfiou(unit,status)
         return
      endif
      extver=0
      hdutype=-1
      call ftghpr(unit,2,simple,bitpix,naxis,naxes,pcount,gcount,extend,status)
      !call ftg2de(unit,0,0.,nmapm,naxes(1),naxes(2),skydld,anyf,status)
      call FTGKYd(unit,'EXPOSURE ',texp,comment,status)
      call ftgkyd(unit,'E1  ',e1,comment,status)
      call ftgkyd(unit,'E2  ',e2,comment,status)
      call ftgkyd(unit,'CRPIX1  ',crpix1,comment,status)
      call ftgkyd(unit,'CRPIX2  ',crpix2,comment,status)
      call ftgkyd(unit,'CRVAL1  ',crval1,comment,status)
      call ftgkyd(unit,'CRVAL2  ',crval2,comment,status)
      call ftgkyd(unit,'CDELT1  ',cdelt1,comment,status)
      call ftgkyd(unit,'CDELT2  ',cdelt2,comment,status)
      call ftgkyd(unit,'CROTA2  ',crota2,comment,status)
      !call ftmnhd(unit,hdutype,'skyerr',extver,status)
      !call ftg2de(unit,0,0.,nmapm,naxes(1),naxes(2),skyerr,anyf,status)

      if(fast .eq. 0) then
         call ftmnhd(unit,hdutype,'skyexp',extver,status)
         call ftg2de(unit,0,0.,nmapm,naxes(1),naxes(2),skyexp,anyf,status)
         skyexp0=skyexp
      else
         skyexp=skyexp0
      endif

      nmapx=naxes(1)
      nmapy=naxes(2)
      call ftclos(unit,status)
      call ftfiou(unit,status)
      if(status.gt.0) then
         return
      endif


      !
      ! Np = (skyexp(n,m)/texp/(e2-e1)) number of open pixels visible from given direction
      ! 
      do n=1,nmapy
         do m=1,nmapx
            
               skyexp(n,m)=skyexp(n,m)/texp/(e2-e1) ! (Np)
            
         enddo
      enddo

      do n=1,nmapx
         do m=1,nmapy
      ! pixel offset
      dist=dsqrt((DBLE(n)-DBLE(ny)/2.0)**2+(DBLE(m)-DBLE(nx)/2.0)**2)
      dist=dist*cdelt1
      image(n,m)=skyexp(n,m)*(1.0+norma0*EXP(-dist/slope0))*(1.0+norma1*EXP(-dist/slope1)) ! Crab fit
         enddo
      enddo
      nx=nmapx
      ny=nmapy

      return
      end