generated from erosita/uds
138 lines
4.2 KiB
Fortran
138 lines
4.2 KiB
Fortran
!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
|
|
|