generated from erosita/uds
236 lines
6.5 KiB
Fortran
236 lines
6.5 KiB
Fortran
!
|
|
! Makes projected GRXE sky model. Model parameters are from Revnivtsev+ 2006
|
|
! The final map can be converted to erg/s/cm2 by multiplying by Lsun/kpc^2
|
|
! Roman Krivonos, 2013
|
|
!
|
|
! ergs-to-photons coefficient is taken from the following model:
|
|
!========================================================================
|
|
!Model atable{polarmodel.fits}<1> + gaussian<2> Source No.: 1 Active/On
|
|
!Model Model Component Parameter Unit Value
|
|
! par comp
|
|
! 1 1 polarcol wdmass 0.600000 frozen
|
|
! 2 1 polarcol norm 1.26483E-09 +/- 2.76500E-11
|
|
! 3 2 gaussian LineE keV 6.70000 frozen
|
|
! 4 2 gaussian Sigma keV 0.190000 frozen
|
|
! 5 2 gaussian norm 1.10106E-02 +/- 1.41770E-03
|
|
!________________________________________________________________________
|
|
! flux 3 20
|
|
! Model Flux 0.22374 photons (2.3232e-09 ergs/cm^2/s) range (3.0000 - 20.000 keV)
|
|
! k= 0.22374/2.3232e-09=9.63068e+07
|
|
|
|
include 'inc/wimage_simple.f'
|
|
include 'inc/deletefile.f'
|
|
include 'inc/printerror.f'
|
|
include 'inc/read_cobe.f'
|
|
|
|
implicit none
|
|
integer nmapm,nmapx,nmapy
|
|
parameter(nmapm=40000)
|
|
integer NROWS
|
|
parameter (NROWS=609)
|
|
real skydld(nmapm,nmapm),skyerr(nmapm,nmapm)
|
|
real skyexp(nmapm,nmapm)
|
|
real sky(nmapm,nmapm)
|
|
|
|
|
|
real*8 crval1,crval2,crpix1,crpix2,cdelt1,cdelt2,crota2
|
|
character*80 comment
|
|
integer unit,status,block,extver,hdutype,simple,bitpix
|
|
integer naxes(3),naxis,pcount,gcount,extend
|
|
logical anyf
|
|
|
|
real*8 ibx1,iby1,ibx2,iby2
|
|
|
|
integer TROWMAX,isel,tnrows,iraw
|
|
parameter (TROWMAX=500000)
|
|
integer sel(TROWMAX)
|
|
real dflux(TROWMAX),dflux_49(TROWMAX),dglat(TROWMAX),dglon(TROWMAX)
|
|
|
|
character filename*140,fileout*140,nullstr*1
|
|
logical exact,anynull
|
|
integer irx,iry,nr
|
|
integer s,i,j,n,m
|
|
real*8 ra,dec,lon,lat,xsrc,ysrc,totalcobe,totalmodel
|
|
real*8 dr,rd,resp2Crab,delta
|
|
character*400 arg
|
|
parameter (dr=0.0174533,rd=57.2958)
|
|
real*8 sumd,sumd_r,sumd_m,koeff,id_l,id_b,b0,texp,sigma,pi,fwhm
|
|
PARAMETER(pi=3.14159)
|
|
real*8 rmin,rmax,rstep,r,rEarth
|
|
real*8 lmin,lmax,bmin,bmax,model,model0,model1,modelrxte
|
|
real*8 dOmega,dV,l,b,x0,y0,z0,theta,phi,zEarth,R0,theta0,phi0
|
|
! parameters for AIT projection
|
|
character*200 outfile
|
|
integer size_sky,sx,sy,k,ir
|
|
parameter (size_sky=2000)
|
|
real ncatched(size_sky,size_sky)
|
|
real*8 xrval,yrval,xinc,yinc,rot,xrpix,yrpix,lcur,bcur,x,y,totFlux,Flux
|
|
real aitmap(nmapm,nmapm),prev,lonr,latr
|
|
|
|
delta=2.0
|
|
sx=int(380.0/delta)
|
|
sy=int(200.0/delta)
|
|
if(sx .gt. nmapm .or. sy .gt. nmapm) then
|
|
print*,sx,sy,nmapm
|
|
return
|
|
endif
|
|
|
|
xinc=delta
|
|
yinc=delta
|
|
rot=0.0
|
|
xrpix=sx/2.0
|
|
yrpix=sy/2.0
|
|
xrval=0.0d0
|
|
yrval=0.0d0
|
|
|
|
|
|
aitmap=0.0
|
|
|
|
zEarth=16.46d0/1000.0d0
|
|
rEarth=8.5d0
|
|
|
|
rmin=1.0
|
|
rmax=22.0
|
|
nr=2600
|
|
rstep=(rmax-rmin)/(nr-1)
|
|
|
|
status=0
|
|
|
|
!lmin=-180.0
|
|
!lmax=180.0
|
|
bmin=-30.0
|
|
bmax=30.0
|
|
|
|
|
|
totFlux=0.0
|
|
sumd=0.0
|
|
sumd_r=0.0
|
|
|
|
outfile='modelrxte_ait_wide.dat'
|
|
call ftgiou(unit,status)
|
|
open(unit,file=outfile)
|
|
|
|
do j=1,nmapm
|
|
do k=1,nmapm
|
|
status=0
|
|
call FTWLDP(DBLE(j),DBLE(k),xrval,yrval,xrpix,yrpix,xinc,yinc,rot,'-AIT',lon,lat,status)
|
|
if(status.eq.0) then
|
|
if(lat.gt.bmax.or.lat.lt.bmin) cycle
|
|
l=lon
|
|
b=lat
|
|
|
|
phi=l*dr
|
|
theta=(90.0+b)*dr
|
|
dOmega=sin(theta)*(delta*dr)**2
|
|
sumd_r=0.0d0
|
|
do ir=0,nr-1
|
|
r=rmin+ir*rstep
|
|
dV=r*r*rstep*dOmega
|
|
x0=r*sin(theta)*cos(phi)-rEarth
|
|
y0=r*sin(theta)*sin(phi)
|
|
z0=r*cos(theta)-zEarth
|
|
R0=dsqrt(x0**2+y0**2+z0**2)
|
|
phi0=datan2(y0,x0)!-pi..pi
|
|
theta0=dacos(z0/R0)!0..pi
|
|
model=modelrxte(x0,y0,z0,R0,theta0,phi0)
|
|
Flux=model*dV/r/r/4/pi
|
|
sumd_r=sumd_r+Flux
|
|
sumd=sumd+model*dV
|
|
totFlux=totFlux+Flux
|
|
enddo
|
|
|
|
aitmap(j,k)=sumd_r
|
|
|
|
write(unit,'(2f9.3,1pe14.4)'),lon,lat,aitmap(j,k)
|
|
else
|
|
call ftgerr(status,comment)
|
|
!print *,comment
|
|
!stop
|
|
cycle
|
|
status=0
|
|
endif
|
|
enddo
|
|
enddo
|
|
|
|
close(unit)
|
|
call ftfiou(unit,status)
|
|
|
|
print*,'Total Flux',totFlux
|
|
print*,'Integral ',sumd
|
|
|
|
! In oder to convert pixel value of the map to flux in this pixel, you need to multiply by Lsun/kpc**2=4.01831e-10
|
|
! Lsun=3.826e+33 erg/s, kpc=3.08568e+21 cm
|
|
|
|
outfile='modelrxte_ait_3to20keV_flux.fits'
|
|
call deletefile(outfile,status)
|
|
call write_image_xy(outfile,aitmap,xrval,yrval,xrpix,yrpix,xinc,yinc,rot,sx,sy,nmapm)
|
|
|
|
200 end
|
|
|
|
|
|
|
|
|
|
real*8 function modelrxte(X,Y,Z,R,theta,phi)
|
|
implicit none
|
|
real*8 dr,rd
|
|
parameter (dr=0.0174533,rd=57.2958)
|
|
real*8, intent(in) :: R,theta,phi
|
|
REAL*8 X,Y,Z
|
|
REAL*8 Xt,Yt,Zt
|
|
REAL*8 X0,Y0,Z0,Rdisk,Rt,Rmax,Rm,BarRmax,Rxy,r0,BarCutoff
|
|
REAL*8 BarTilt,Zdisk,rhobulge0,rhobulge,rhodisk,rhodisk0
|
|
parameter (BarTilt=(180d0-29.0)*dr,BarRmax=2.4)
|
|
parameter (X0=3.4,Y0=1.2,Z0=1.12,Rdisk=2.5,Zdisk=0.13,Rm=3.0)
|
|
real*8 Lsun
|
|
parameter(Lsun = 3.826e+33)
|
|
|
|
modelrxte=0.0
|
|
|
|
|
|
rhobulge0=1.0
|
|
! Mikej bulge luminosity 3.9+/-0.5 x 10^37 erg/s = 1.02e4 Lsun
|
|
! bulge volume without BarCutoff = 42.599
|
|
rhobulge0=1.02e4/42.599
|
|
|
|
Rmax=10.0
|
|
|
|
rhodisk0=1.0
|
|
! disk volume = 4.82471
|
|
! Mikej disk lumin ~10e37 = 26137 Lsun
|
|
rhodisk0=2.6137e4/4.82471
|
|
|
|
Xt=X*dcos(BarTilt)-Y*dsin(BarTilt)
|
|
Yt=X*dsin(BarTilt)+Y*dcos(BarTilt)
|
|
Zt=Z
|
|
|
|
r0=0.5
|
|
Rxy=dsqrt( (Xt)**2 + (Yt)**2 )
|
|
BarCutoff=1 .0
|
|
if(Rxy.gt.BarRmax) then
|
|
BarCutoff=dexp(-Rxy**2/2/r0**2)
|
|
endif
|
|
|
|
|
|
Rt=dsqrt((Xt/X0)**2 + (Yt/Y0)**2 + (Zt/Z0)**2)
|
|
|
|
! Avoid singularity in GC, say close than 1pc
|
|
if(Rt.lt.0.05) Rt=0.05
|
|
|
|
rhobulge = rhobulge0*Rt**(-1.8)*exp(-Rt**3)*BarCutoff
|
|
|
|
rhodisk = 0.0
|
|
if(R.lt.Rmax) then
|
|
rhodisk=rhodisk0*dexp(-(Rm/R)**3 -R/Rdisk -dabs(Z)/Zdisk)
|
|
endif
|
|
|
|
modelrxte=rhodisk+rhobulge
|
|
|
|
|
|
!check, volume should be 4/3*pi = 4.18879
|
|
!modelrxte=0.0
|
|
!if(R.lt.15.0) modelrxte=1.0
|
|
|
|
return
|
|
end function modelrxte
|
|
|