diff --git a/data/allsky/Makefile b/data/allsky/Makefile new file mode 100644 index 0000000..6ffb680 --- /dev/null +++ b/data/allsky/Makefile @@ -0,0 +1,10 @@ +LIBS=-L ${HEADAS}/lib -lcfitsio + +FFLAGS= -O -ffree-form + +all:modelrxte + +modelrxte: modelrxte.f + ${FC} ${FFLAGS} -o modelrxte modelrxte.f ${LIBS} + + diff --git a/data/allsky/inc/constants.inc b/data/allsky/inc/constants.inc new file mode 100644 index 0000000..7b64314 --- /dev/null +++ b/data/allsky/inc/constants.inc @@ -0,0 +1,60 @@ +c physical constants +c MG; Mon Nov 11 16:28:45 MSK 1996 + + +c------ pi + real*8 pi,pi2,rg,gr + parameter (pi=3.14159265358979312d00,pi2=pi+pi,rg=180d0/pi,gr=1d0/rg) + +c------ physical constants + real*8 sigthomson + parameter(sigthomson=6.65246d-25) + real*8 clight + parameter(clight=2.99792458d10) + real*8 mproton,melectron + parameter(mproton=1.6726231d-24,melectron=9.1093897d-28) + real*8 echarge + parameter(echarge=4.8032068d-10) + real*8 sigsb + parameter(sigsb=5.67051d-5) + real*8 kboltzmann + parameter(kboltzmann=1.380658d-16) + real*8 ggrav + parameter(ggrav=6.67259d-8) + real*8 hplanck + parameter(hplanck=6.6260755d-27) + real*8 msun + parameter (msun=1.9891d33) + real*8 rsun + parameter (rsun=6.9599d10) + real*8 lsun + parameter (lsun=3.845d33) + + + + + + + +c------ units conversions +c...... energy + real*8 kev2erg + parameter(kev2erg=1.602177d-9) + real*8 kev2hz + parameter(kev2hz=2.417965d17) + real*8 kev2kelvin + parameter(kev2kelvin=1.16048d7) + real*8 jy2erg + parameter(jy2erg=1d-23) +c...... flux +c flux[erg/s/cm^2/Hz] = ergcms_anstr2ergcms_hz_0 * lambda[Anstr]^2 * flux[erg/s/cm^2/A] + real*8 ergcms_anstr2ergcms_hz_0 + parameter(ergcms_anstr2ergcms_hz_0=3.33d-19) +c...... lenth + real*8 pc2cm,kpc2cm,mpc2cm + parameter(pc2cm=3.085678d18,kpc2cm=1e3*pc2cm,mpc2cm=1e6*pc2cm) +c...... time + real*8 yr2sec + parameter(yr2sec=3.15576d7) + + diff --git a/data/allsky/inc/crab.f b/data/allsky/inc/crab.f new file mode 100644 index 0000000..e316430 --- /dev/null +++ b/data/allsky/inc/crab.f @@ -0,0 +1,21 @@ +!c========================================================================== + subroutine printerror(status) +!c-------------------------------------------------------------------------- +!c print error message (for fitsio routines) +!c-------------------------------------------------------------------------- + implicit none + integer status + character errtext*30,errmessage*80 + + if(status.le.0) return + call ftgerr(status,errtext) + print *,'FITSIO Error Status =',status,': ',errtext + call ftgmsg(errmessage) + do while (errmessage .ne. ' ') + print *,errmessage + call ftgmsg(errmessage) + end do +!c if(status.ne.0) stop + + return + end diff --git a/data/allsky/inc/crabmodel.f b/data/allsky/inc/crabmodel.f new file mode 100644 index 0000000..5fe5b3d --- /dev/null +++ b/data/allsky/inc/crabmodel.f @@ -0,0 +1,21 @@ +!c========================================================================== + subroutine crabmodel(status) +!c-------------------------------------------------------------------------- +!c print error message (for fitsio routines) +!c-------------------------------------------------------------------------- + implicit none + integer status + character errtext*30,errmessage*80 + + if(status.le.0) return + call ftgerr(status,errtext) + print *,'FITSIO Error Status =',status,': ',errtext + call ftgmsg(errmessage) + do while (errmessage .ne. ' ') + print *,errmessage + call ftgmsg(errmessage) + end do +!c if(status.ne.0) stop + + return + end diff --git a/data/allsky/inc/deletefile.f b/data/allsky/inc/deletefile.f new file mode 100644 index 0000000..370c410 --- /dev/null +++ b/data/allsky/inc/deletefile.f @@ -0,0 +1,20 @@ +!c========================================================================= + subroutine deletefile(filename,status) +!c------------------------------------------------------------------------- +!c delete fits file +!c------------------------------------------------------------------------- + implicit none + integer status + character*(*) filename + + integer access,i,unlink + + + + + if(access(filename,'w').eq.0) then + i=unlink(filename) + endif + + return + end diff --git a/data/allsky/inc/ibis_eresp.f b/data/allsky/inc/ibis_eresp.f new file mode 100644 index 0000000..c9b5a04 --- /dev/null +++ b/data/allsky/inc/ibis_eresp.f @@ -0,0 +1,137 @@ +!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 + diff --git a/data/allsky/inc/printerror.f b/data/allsky/inc/printerror.f new file mode 100644 index 0000000..e316430 --- /dev/null +++ b/data/allsky/inc/printerror.f @@ -0,0 +1,21 @@ +!c========================================================================== + subroutine printerror(status) +!c-------------------------------------------------------------------------- +!c print error message (for fitsio routines) +!c-------------------------------------------------------------------------- + implicit none + integer status + character errtext*30,errmessage*80 + + if(status.le.0) return + call ftgerr(status,errtext) + print *,'FITSIO Error Status =',status,': ',errtext + call ftgmsg(errmessage) + do while (errmessage .ne. ' ') + print *,errmessage + call ftgmsg(errmessage) + end do +!c if(status.ne.0) stop + + return + end diff --git a/data/allsky/inc/read_co.f b/data/allsky/inc/read_co.f new file mode 100644 index 0000000..eb7105e --- /dev/null +++ b/data/allsky/inc/read_co.f @@ -0,0 +1,110 @@ +!c=============================================================== + subroutine read_co(dglon,dglat,dflux,tnrows) +!c-------------------------------------------------------------- +!c write image to FITS file +!c +!c Created: Mon May 16 14:21:02 EDT 1994 +!c-------------------------------------------------------------- + implicit none + character filename*40 + integer nxm,nx,ny,erange,n,m + logical anyf + logical exact,anynull + + integer lnblnk + 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) + integer nhdu,tnrows,tncols + integer ixc,iyc,ival + real skyexp(nmapm,nmapm) + real skydld(nmapm,nmapm) + real skyerr(nmapm,nmapm) + real image(nmapm,nmapm) + + ! cobe + integer TROWMAX,ipixno,irezid,istddev,iglat,iglon,iraw,isum,isel + parameter (TROWMAX=500000) + integer dpixno(TROWMAX),idx,sel(TROWMAX) + real dflux(TROWMAX),dglat(TROWMAX),dglon(TROWMAX) + real*8 dl,db + + real dflux_49(TROWMAX) + integer*4 mask(TROWMAX) + + + + call ftgiou(unit,status) + call ftopen(unit,'/afs/mpa/home/krivonos/results/latscan/CO/comap_list.fits',0,block,status) + if(status.ne.0) then + print*,'Bad file, skip it' + call ftclos(unit,status) + call ftfiou(unit,status) + status=0 + call exit(0) + endif + call FTGHDN(unit,nhdu) + call ftmahd(unit,2,hdutype,status) + call FTGHDN(unit,nhdu) + call FTGNRW(unit,tnrows, status) + call FTGNCL(unit,tncols, status) + + if(tnrows.gt.trowmax) then + print*,'Too many points',tnrows,trowmax + call ftclos(unit,status) + call ftfiou(unit,status) + stop + endif + exact=.false. + call ftgcno(unit,exact,'XCOORD',ixc,status) + call ftgcno(unit,exact,'YCOORD',iyc,status) + call ftgcno(unit,exact,'VALUE',ival,status) + + do iraw=1,tnrows + call ftgcve(unit,ixc,iraw,1,1,0,dglon(iraw),anynull,status) + call ftgcve(unit,iyc,iraw,1,1,0,dglat(iraw),anynull,status) + call ftgcve(unit,ival,iraw,1,1,0.,dflux(iraw),anynull,status) + end do + + call ftclos(unit,status) + call ftfiou(unit,status) + if(status.gt.0) then + print*,'Too many errors',tnrows,trowmax + call ftclos(unit,status) + call ftfiou(unit,status) + call exit(0) + endif + + + do iraw=1,tnrows + if (isnan(dflux(iraw))) then + !if(abs(Nint(dflux(iraw))) .eq. 32768) then + print*,dglon(iraw),dglat(iraw),dflux(iraw) + dflux(iraw)=0.0 + endif + enddo + + ! modify coordinates + do iraw=1,tnrows + if(dglon(iraw) .gt. 180.0) then + dglon(iraw)=dglon(iraw)-360.0 + endif + enddo + + + !RAW MAP + !do iraw=1,tnrows + ! print*,dglon(iraw),dglat(iraw),dflux(iraw) + !enddo + !stop + + + + return + end + diff --git a/data/allsky/inc/read_cobe.f b/data/allsky/inc/read_cobe.f new file mode 100644 index 0000000..1f4bb6d --- /dev/null +++ b/data/allsky/inc/read_cobe.f @@ -0,0 +1,219 @@ +!c=============================================================== + subroutine read_cobe(dglon,dglat,dflux_49,tnrows) +!c-------------------------------------------------------------- +!c write image to FITS file +!c +!c Created: Mon May 16 14:21:02 EDT 1994 +!c-------------------------------------------------------------- + implicit none + character filename*40 + integer nxm,nx,ny,erange,n,m + logical anyf + logical exact,anynull + + integer lnblnk + 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) + integer nhdu,tnrows,tncols + + real skyexp(nmapm,nmapm) + real skydld(nmapm,nmapm) + real skyerr(nmapm,nmapm) + real image(nmapm,nmapm) + + ! cobe + integer TROWMAX,ipixno,irezid,istddev,iglat,iglon,iraw,isum,isel + parameter (TROWMAX=500000) + integer dpixno(TROWMAX),idx,sel(TROWMAX) + real dflux(TROWMAX),dglat(TROWMAX),dglon(TROWMAX) + real*8 dl,db + + real dflux_49(TROWMAX),dflux_49_orig(TROWMAX),dflux_12(TROWMAX) + real dratio(TROWMAX) + integer*4 mask(TROWMAX) + +! read DIRBE flux 4.9 microns (Resid) +! The sum of the selected column is 90891.360 +! The mean of the selected column is 0.23114868 +! The standard deviation of the selected column is 0.76384730 +! The minimum of selected column is 3.64903696E-02 +! The maximum of selected column is 110.34200 +! The number of points used in calculation is 393216 +! +! Expected value 0.0823 +! + + call ftgiou(unit,status) + call ftopen(unit,'/afs/mpa/project/integral/results/latscan/COBE/DIRBE_BAND04_ZSMA.FITS',0,block,status) + !call ftopen(unit,'/afs/mpa/project/integral/results/latscan/COBE/DIRBE_BAND3A_ZSMA.FITS',0,block,status) + if(status.ne.0) then + print*,'Bad file, skip it' + call ftclos(unit,status) + call ftfiou(unit,status) + status=0 + call exit(0) + endif + call FTGHDN(unit,nhdu) + call ftmahd(unit,2,hdutype,status) + call FTGHDN(unit,nhdu) + call FTGNRW(unit,tnrows, status) + call FTGNCL(unit,tncols, status) + + if(tnrows.gt.trowmax) then + print*,'Too many points',tnrows,trowmax + call ftclos(unit,status) + call ftfiou(unit,status) + stop + endif + exact=.false. + call ftgcno(unit,exact,'pixel_no',ipixno,status) + call ftgcno(unit,exact,'resid',irezid,status) + + do iraw=1,tnrows + call ftgcvj(unit,ipixno,iraw,1,1,0,dpixno(iraw),anynull,status) + call ftgcve(unit,irezid,iraw,1,1,0.,dflux_49(iraw),anynull,status) + end do + + call ftclos(unit,status) + call ftfiou(unit,status) + if(status.gt.0) then + call exit(0) + endif + +! +! Read DIRBE flux 1.25 microns +! +! The sum of the selected column is 194936.15 +! The mean of the selected column is 0.49574827 +! The standard deviation of the selected column is 1.3725306 +! The minimum of selected column is 2.56427769E-02 +! The maximum of selected column is 215.81805 +! The number of points used in calculation is 393216 +! +! Expected value: 0.0775 +! + call ftgiou(unit,status) + call ftopen(unit,'/afs/mpa/project/integral/results/latscan/COBE/DIRBE_BAND1A_ZSMA.FITS',0,block,status) + if(status.ne.0) then + print*,'Bad file, skip it' + call ftclos(unit,status) + call ftfiou(unit,status) + status=0 + call exit(0) + endif + call FTGHDN(unit,nhdu) + call ftmahd(unit,2,hdutype,status) + call FTGHDN(unit,nhdu) + call FTGNRW(unit,tnrows, status) + call FTGNCL(unit,tncols, status) + + if(tnrows.gt.trowmax) then + print*,'Too many points',tnrows,trowmax + call ftclos(unit,status) + call ftfiou(unit,status) + stop + endif + exact=.false. + call ftgcno(unit,exact,'pixel_no',ipixno,status) + call ftgcno(unit,exact,'resid',irezid,status) + + do iraw=1,tnrows + call ftgcvj(unit,ipixno,iraw,1,1,0,dpixno(iraw),anynull,status) + call ftgcve(unit,irezid,iraw,1,1,0.,dflux_12(iraw),anynull,status) + end do + + call ftclos(unit,status) + call ftfiou(unit,status) + if(status.gt.0) then + call exit(0) + endif + + ! read DIRBE coordinates + call ftgiou(unit,status) + call ftopen(unit,'/afs/mpa/project/integral/results/latscan/COBE/DIRBE_SKYMAP_INFO.FITS',0,block,status) + if(status.ne.0) then + print*,'Bad file, skip it' + call ftclos(unit,status) + call ftfiou(unit,status) + status=0 + call exit(0) + endif + call FTGHDN(unit,nhdu) + call ftmahd(unit,2,hdutype,status) + call FTGHDN(unit,nhdu) + call FTGNRW(unit,tnrows, status) + call FTGNCL(unit,tncols, status) + + if(tnrows.gt.trowmax) then + print*,'Too many points',tnrows,trowmax + call ftclos(unit,status) + call ftfiou(unit,status) + stop + endif + exact=.false. + call ftgcno(unit,exact,'QSPIXEL',ipixno,status) + call ftgcno(unit,exact,'GLON-CSC',iglon,status) + call ftgcno(unit,exact,'GLAT-CSC',iglat,status) + + do iraw=1,tnrows + call ftgcvj(unit,ipixno,iraw,1,1,0,dpixno(iraw),anynull,status) + call ftgcve(unit,iglon,iraw,1,1,0.,dglon(iraw),anynull,status) + call ftgcve(unit,iglat,iraw,1,1,0.,dglat(iraw),anynull,status) + end do + + call ftclos(unit,status) + call ftfiou(unit,status) + if(status.gt.0) then + call exit(0) + endif + + ! modify coordinates + do iraw=1,tnrows + if(dglon(iraw) .gt. 180.0) then + dglon(iraw)=dglon(iraw)-360.0 + endif + enddo + + ! substruct background + do iraw=1,tnrows + dflux_12(iraw)=dflux_12(iraw)-0.07 + dflux_49(iraw)=dflux_49(iraw)-0.08 + enddo + + ! correction for absorbtion + do iraw=1,tnrows + dratio(iraw)=0.0 + if(dflux_49(iraw).gt.0.5.and.dflux_49(iraw).lt.50.0.and.dflux_12(iraw).gt. 0.8 .and. dflux_12(iraw).lt.20.0) then + dratio(iraw)=dflux_12(iraw)/dflux_49(iraw)/3.5 + else + dflux_49(iraw)=0.0 + dflux_12(iraw)=0.0 + dratio(iraw)=0.0 + endif + enddo + + do iraw=1,tnrows + if(dratio(iraw) .gt. 0.0 .and. dratio(iraw) .lt. 1.0) then + dratio(iraw)=(dratio(iraw))**(-0.25) + dflux_49(iraw)=dratio(iraw)*dflux_49(iraw) + else + dflux_49(iraw)=0.0 + endif + enddo + + do iraw=1,tnrows + if(abs(dglat(iraw)).gt.7.0) then + dflux_49(iraw)=0.0 + endif + enddo + + + return + end + diff --git a/data/allsky/inc/read_modelrxte.f b/data/allsky/inc/read_modelrxte.f new file mode 100644 index 0000000..bc9da85 --- /dev/null +++ b/data/allsky/inc/read_modelrxte.f @@ -0,0 +1,104 @@ +!c=============================================================== + subroutine read_modelrxte(dglon,dglat,dflux,tnrows) +!c-------------------------------------------------------------- +!c write image to FITS file +!c +!c Created: Mon May 16 14:21:02 EDT 1994 +!c-------------------------------------------------------------- +! see details /afs/mpa/project/integral/results/latscan/f90/modelrxte.f + implicit none + character filename*40 + integer nxm,nx,ny,erange,n,m + logical anyf + logical exact,anynull + integer lnblnk + 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,fn + integer nmapm,nmapx,nmapy + parameter(nmapm=2000) + integer nhdu,tnrows,tncols + integer ixc,iyc,ival + real skyexp(nmapm,nmapm) + real skydld(nmapm,nmapm) + real skyerr(nmapm,nmapm) + real image(nmapm,nmapm) + integer lfblnk + ! cobe + integer TROWMAX,ipixno,irezid,istddev,iglat,iglon,iraw,isum,isel + parameter (TROWMAX=3000000) + integer dpixno(TROWMAX),idx,sel(TROWMAX) + real dflux(TROWMAX),dglat(TROWMAX),dglon(TROWMAX) + real*8 dl,db + + !real dflux_49(TROWMAX) + integer*4 mask(TROWMAX) + + status=0 + !fn='/afs/mpa/project/integral/results/latscan/f90/modelrxte_hires.bintab.fits' + fn='/afs/mpa/project/integral/results/latscan/f90/modelrxte_ait_wide.bintab.fits' + + print*,'Read file',fn(:lfblnk(fn)) + call ftgiou(unit,status) + call ftopen(unit,fn,0,block,status) + if(status.ne.0) then + print*,'Bad file, skip it' + call ftclos(unit,status) + call ftfiou(unit,status) + status=0 + call exit(0) + endif + call FTGHDN(unit,nhdu) + call ftmahd(unit,2,hdutype,status) + call FTGHDN(unit,nhdu) + call FTGNRW(unit,tnrows, status) + call FTGNCL(unit,tncols, status) + + if(tnrows.gt.trowmax) then + print*,'Too many points',tnrows,trowmax + call ftclos(unit,status) + call ftfiou(unit,status) + stop + endif + exact=.false. + call ftgcno(unit,exact,'GLON',ixc,status) + call ftgcno(unit,exact,'GLAT',iyc,status) + call ftgcno(unit,exact,'VALUE',ival,status) + + do iraw=1,tnrows + call ftgcve(unit,ixc,iraw,1,1,0,dglon(iraw),anynull,status) + call ftgcve(unit,iyc,iraw,1,1,0,dglat(iraw),anynull,status) + call ftgcve(unit,ival,iraw,1,1,0.,dflux(iraw),anynull,status) + end do + + call ftclos(unit,status) + call ftfiou(unit,status) + if(status.gt.0) then + print*,'Too many errors',tnrows,trowmax + call ftclos(unit,status) + call ftfiou(unit,status) + call exit(0) + endif + + + do iraw=1,tnrows + if (isnan(dflux(iraw))) then + !if(abs(Nint(dflux(iraw))) .eq. 32768) then + print*,dglon(iraw),dglat(iraw),dflux(iraw) + dflux(iraw)=0.0 + endif + enddo + + ! modify coordinates + do iraw=1,tnrows + if(dglon(iraw) .gt. 180.0) then + dglon(iraw)=dglon(iraw)-360.0 + endif + enddo + + return + end + diff --git a/data/allsky/inc/wimage_simple.f b/data/allsky/inc/wimage_simple.f new file mode 100644 index 0000000..9076555 --- /dev/null +++ b/data/allsky/inc/wimage_simple.f @@ -0,0 +1,279 @@ +!c=============================================================== + subroutine wimage_simple_i4(filename,image,nxm,nx,ny) +!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 + integer image(nxm,*) + + integer lnblnk + integer status,bitpix,naxis,pcount,gcount + integer naxes(2) + integer unit + + status=0 + print 10,filename(:lnblnk(filename)) + 10 format(/'Write FITS image to file ',a) + + call deletefile(filename,status) + call ftgiou(unit,status) + + bitpix=32 + naxis=2 + pcount=0 + gcount=1 + naxes(1)=nx + naxes(2)=ny + + call ftinit(unit,filename,0,status) + call ftphpr(unit,.true.,bitpix,naxis,naxes,pcount,gcount,.false.,status) + call ftpdef(unit,bitpix,naxis,naxes,pcount,gcount,status) + call ftp2dj(unit,0,nxm,naxes(1),naxes(2),image,status) + call ftclos(unit,status) + call ftfiou(unit,status) + if(status.gt.0) call printerror(status) + + + return + end +!c=============================================================== + subroutine wimage_simple_r4(filename,image,nxm,nx,ny) +!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 + real image(nxm,*) + + integer lnblnk + integer status,bitpix,naxis,pcount,gcount + integer naxes(2) + integer unit + + status=0 + print 10,filename(:lnblnk(filename)) + 10 format(/'Write FITS image to file ',a) + + call deletefile(filename,status) + call ftgiou(unit,status) + + bitpix=-32 + naxis=2 + pcount=0 + gcount=1 + naxes(1)=nx + naxes(2)=ny + + call ftinit(unit,filename,0,status) + call ftphpr(unit,.true.,bitpix,naxis,naxes,pcount,gcount,.false.,status) + call ftpdef(unit,bitpix,naxis,naxes,pcount,gcount,status) + call ftp2de(unit,0,nxm,naxes(1),naxes(2),image,status) + call ftclos(unit,status) + call ftfiou(unit,status) + if(status.gt.0) call printerror(status) + + + return + end +!c=============================================================== + subroutine wimage_simple_r8(filename,image,nxm,nx,ny) +!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 + real*8 image(nxm,*) + + integer lnblnk + integer status,bitpix,naxis,pcount,gcount + integer naxes(2) + integer unit + + status=0 + print 10,filename(:lnblnk(filename)) + 10 format(/'Write FITS image to file ',a) + + call deletefile(filename,status) + call ftgiou(unit,status) + + bitpix=-64 + naxis=2 + pcount=0 + gcount=1 + naxes(1)=nx + naxes(2)=ny + + call ftinit(unit,filename,0,status) + call ftphpr(unit,.true.,bitpix,naxis,naxes,pcount,gcount,.false.,status) + call ftpdef(unit,bitpix,naxis,naxes,pcount,gcount,status) + call ftp2dd(unit,0,nxm,naxes(1),naxes(2),image,status) + call ftclos(unit,status) + call ftfiou(unit,status) + if(status.gt.0) call printerror(status) + + + return + end +!c============================================================================= + subroutine read_image_r4(fn,image,nxm,mapx,mapy) +!c----------------------------------------------------------------------------- +!c read image (real*4) +!c----------------------------------------------------------------------------- + implicit none + character*(*) fn + integer nxm,mapx,mapy + real image(nxm,*) + + character*80 comment + integer unit,status,block,hdutype,extver,naxis,naxes(4),bitpix + integer pcount,gcount + logical anyf,exact,simple,extend + + status=0 + call ftgiou(unit,status) + call ftopen(unit,fn,0,block,status) + + call ftghpr(unit,2,simple,bitpix,naxis,naxes,pcount,gcount,extend,status) + !print*,' image nx,ny:',naxes(1),naxes(2) + call ftg2de(unit,0,0.,nxm,naxes(1),naxes(2),image,anyf,status) + call ftclos(unit,status) + call ftfiou(unit,status) + if(status.gt.0) then + call printerror(status) + endif + + mapx=naxes(1) + mapy=naxes(2) + + return + end + + + subroutine write_image(outfile,flux,xrval,yrval,xrpix,yrpix,xinc,yinc,rot,size_sky,nxm) + implicit none + integer nxm + integer size_sky,i,j + real*8 xrval,yrval,xrpix,yrpix,xinc,yinc,rot + real flux(nxm,*) + + + integer ounit,status,lnblnk + integer bitpix,naxis,naxes(2),pcount,gcount + character*200 outfile,errmsg + + bitpix=-32 + naxis=2 + pcount=0 + gcount=1 + naxes(1)=size_sky + naxes(2)=size_sky + + + + status=0 + call ftgiou(ounit,status) + call ftinit(ounit,outfile,0,status) + + if(status.eq.103) then + print*,' Can not create file. File already exist' + print *,'' + stop + endif + call ftphpr(ounit,.true.,bitpix,naxis,naxes,pcount,gcount,.true.,status) + call ftpdef(ounit,bitpix,naxis,naxes,pcount,gcount,status) + + call ftpkys(ounit,'CTYPE1','GLON-AIT','X-axis type',status) + call ftpkys(ounit,'CTYPE2','GLAT-AIT','Y-axis type',status) + call ftpkyd(ounit,'CRVAL1',xrval,10,'Reference pixel value',status) + call ftpkyd(ounit,'CRVAL2',yrval,10,'Reference pixel value',status) + call ftpkyd(ounit,'CRPIX1',xrpix,10,'Reference pixel',status) + call ftpkyd(ounit,'CRPIX2',yrpix,10,'Reference pixel',status) + call ftpkyd(ounit,'CDELT1',xinc,10,'Degrees/pixel',status) + call ftpkyd(ounit,'CDELT2',yinc,10,'Degrees/pixel',status) + call ftpkyd(ounit,'CROTA1',rot,10,'',status) + call ftpdat(ounit,status) + call ftphis(ounit,'Created by R.Krivonos, krivonos@iki.rssi.ru',status) + + call ftp2de(ounit,0,nxm,size_sky,size_sky,flux,status) + if (status.ne.0) then + call ftgerr(status,errmsg) + print *,errmsg + endif + call ftpcks(ounit,status) + call ftclos(ounit,status) + call ftfiou(ounit,status) + if(status.eq.0) then + print *,'File is written: ',outfile(1:lnblnk(outfile)) + endif + + end + + subroutine write_image_xy(outfile,flux,xrval,yrval,xrpix,yrpix,xinc,yinc,rot,size_sky_x,size_sky_y,nxm) + implicit none + integer nxm + integer size_sky_x,size_sky_y,i,j + real*8 xrval,yrval,xrpix,yrpix,xinc,yinc,rot + real flux(nxm,*) + + + integer ounit,status,lnblnk + integer bitpix,naxis,naxes(2),pcount,gcount + character*200 outfile,errmsg + + bitpix=-32 + naxis=2 + pcount=0 + gcount=1 + naxes(1)=size_sky_x + naxes(2)=size_sky_y + + + + status=0 + call ftgiou(ounit,status) + call ftinit(ounit,outfile,0,status) + + if(status.eq.103) then + print*,' Can not create file. File already exist' + print *,'' + stop + endif + call ftphpr(ounit,.true.,bitpix,naxis,naxes,pcount,gcount,.true.,status) + call ftpdef(ounit,bitpix,naxis,naxes,pcount,gcount,status) + + call ftpkys(ounit,'CTYPE1','GLON-AIT','X-axis type',status) + call ftpkys(ounit,'CTYPE2','GLAT-AIT','Y-axis type',status) + call ftpkyd(ounit,'CRVAL1',xrval,10,'Reference pixel value',status) + call ftpkyd(ounit,'CRVAL2',yrval,10,'Reference pixel value',status) + call ftpkyd(ounit,'CRPIX1',xrpix,10,'Reference pixel',status) + call ftpkyd(ounit,'CRPIX2',yrpix,10,'Reference pixel',status) + call ftpkyd(ounit,'CDELT1',xinc,10,'Degrees/pixel',status) + call ftpkyd(ounit,'CDELT2',yinc,10,'Degrees/pixel',status) + call ftpkyd(ounit,'CROTA1',rot,10,'',status) + call ftpdat(ounit,status) + call ftphis(ounit,'Created by R.Krivonos',status) + + call ftp2de(ounit,0,nxm,size_sky_x,size_sky_y,flux,status) + if (status.ne.0) then + call ftgerr(status,errmsg) + print *,errmsg + endif + call ftpcks(ounit,status) + call ftclos(ounit,status) + call ftfiou(ounit,status) + if(status.eq.0) then + print *,'File is written: ',outfile(1:lnblnk(outfile)) + endif + + end + diff --git a/data/allsky/modelrxte.f b/data/allsky/modelrxte.f new file mode 100644 index 0000000..b5f9b47 --- /dev/null +++ b/data/allsky/modelrxte.f @@ -0,0 +1,235 @@ +! +! 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 + diff --git a/data/allsky/modelrxte_ait_3to20keV_flux_004deg.fits b/data/allsky/modelrxte_ait_3to20keV_flux_004deg.fits new file mode 100644 index 0000000..4a69f92 Binary files /dev/null and b/data/allsky/modelrxte_ait_3to20keV_flux_004deg.fits differ diff --git a/data/allsky/modelrxte_ait_3to20keV_flux_006deg.fits b/data/allsky/modelrxte_ait_3to20keV_flux_006deg.fits new file mode 100644 index 0000000..1145be7 Binary files /dev/null and b/data/allsky/modelrxte_ait_3to20keV_flux_006deg.fits differ diff --git a/data/allsky/modelrxte_ait_3to20keV_flux_01deg.fits b/data/allsky/modelrxte_ait_3to20keV_flux_01deg.fits new file mode 100644 index 0000000..cca8b6c Binary files /dev/null and b/data/allsky/modelrxte_ait_3to20keV_flux_01deg.fits differ diff --git a/data/allsky/modelrxte_ait_3to20keV_flux_02deg.fits b/data/allsky/modelrxte_ait_3to20keV_flux_02deg.fits new file mode 100644 index 0000000..e43297f Binary files /dev/null and b/data/allsky/modelrxte_ait_3to20keV_flux_02deg.fits differ diff --git a/data/allsky/modelrxte_ait_3to20keV_flux_03deg.fits b/data/allsky/modelrxte_ait_3to20keV_flux_03deg.fits new file mode 100644 index 0000000..d561377 Binary files /dev/null and b/data/allsky/modelrxte_ait_3to20keV_flux_03deg.fits differ diff --git a/data/allsky/modelrxte_ait_3to20keV_flux_04deg.fits b/data/allsky/modelrxte_ait_3to20keV_flux_04deg.fits new file mode 100644 index 0000000..8585793 Binary files /dev/null and b/data/allsky/modelrxte_ait_3to20keV_flux_04deg.fits differ diff --git a/data/modelrxte_ait_3to20keV_flux_1deg.fits b/data/allsky/modelrxte_ait_3to20keV_flux_1deg.fits similarity index 100% rename from data/modelrxte_ait_3to20keV_flux_1deg.fits rename to data/allsky/modelrxte_ait_3to20keV_flux_1deg.fits diff --git a/data/modelrxte_ait_3to20keV_flux_2deg.fits b/data/allsky/modelrxte_ait_3to20keV_flux_2deg.fits similarity index 100% rename from data/modelrxte_ait_3to20keV_flux_2deg.fits rename to data/allsky/modelrxte_ait_3to20keV_flux_2deg.fits diff --git a/data/xspec/LON+20_bremms.xcm b/data/xspec/LON+20_bremms.xcm new file mode 100644 index 0000000..45251a8 --- /dev/null +++ b/data/xspec/LON+20_bremms.xcm @@ -0,0 +1,18 @@ +method leven 10 0.01 +abund angr +xsect vern +cosmo 70 0 0.73 +xset delta 0.01 +systematic 0 +model cflux*bremss + cflux*powerlaw + 20 -0.1 0 0 1e+06 1e+06 + 60 -0.1 0 0 1e+06 1e+06 + -9.0642 0.01 -100 -100 100 100 + 20.3536 1 0.0001 0.0001 100 200 + 1e-05 -1 0 0 1e+20 1e+24 + 80 -0.1 0 0 1e+06 1e+06 + 200 -0.1 0 0 1e+06 1e+06 + -11.2811 0.01 -100 -100 100 100 + -1.5 -1 -3 -2 9 10 + 1e-05 -1 0 0 1e+20 1e+24 +bayes off diff --git a/scripts/02_grxe_resid.py b/scripts/02_grxe_resid.py index 0c50b84..046ade7 100755 --- a/scripts/02_grxe_resid.py +++ b/scripts/02_grxe_resid.py @@ -60,7 +60,7 @@ mjd = data.field('mjd') clean = data.field('clean') phase = data.field('phase') - +obsid0=[] rev0=[] phase0=[] clean0=[] @@ -133,6 +133,7 @@ for i, row in df.iterrows(): r1 = bgdmodel[orbit]['r1'] # nearest left orbit used for calibration r2 = bgdmodel[orbit]['r2'] # nearest right orbit used for calibration + obsid0.append(obsid) c0.append(c) base0.append(abs(orbit - int(np.min([r1,r2])))) clean0.append(clean[i]) @@ -160,7 +161,7 @@ indices = sorted( ) coldefs = fits.ColDefs([ - #fits.Column(name='OBSID', format='11A', array=[obs_id[index] for index in indices]), + fits.Column(name='OBSID', format='11A', array=[obsid0[index] for index in indices]), #fits.Column(name='RA', format='D', unit='deg', array=[ra[index] for index in indices]), #fits.Column(name='DEC', format='D', unit='deg', array=[dec[index] for index in indices]), fits.Column(name='LON', format='D', unit='deg', array=[lon0[index] for index in indices]), diff --git a/scripts/03_grxe_map.py b/scripts/03_grxe_map.py index 30ee53f..13bd62e 100755 --- a/scripts/03_grxe_map.py +++ b/scripts/03_grxe_map.py @@ -56,7 +56,7 @@ minmax=[-300,800] sigma=3 maxiters=10 -modelrxte="modelrxte_ait_3to20keV_flux_2deg.fits" +modelrxte="allsky/modelrxte_ait_3to20keV_flux_2deg.fits" hdulist = fits.open(datadir+modelrxte) w = wcs.WCS(hdulist[0].header) smap =hdulist[0].data @@ -87,6 +87,10 @@ sign_map = np.array([[0.0 for i in range(sx)] for j in range(sy)]) sem_map = np.array([[0.0 for i in range(sx)] for j in range(sy)]) cnt_map = np.array([[0 for i in range(sx)] for j in range(sy)]) +obsid_map = {} +grxe_map = {} +grxe_err_map = {} + for i in range(sx): for j in range(sy): world = w.wcs_pix2world([(i+1,j+1)], 1) @@ -94,6 +98,9 @@ for i in range(sx): lat = world[0][1] if(np.isnan(lon) or np.isnan(lat)): continue + + + ds9i=i+1 ds9j=j+1 df0 = df.query('DS9X == {} & DS9Y == {}'.format(ds9i,ds9j)) @@ -104,8 +111,17 @@ for i in range(sx): #print("***",i+1,j+1,lon,lat,smap[j][i]) #for i0,row0 in df0.iterrows(): # print(row0['LON'],row0['LAT'],row0['GRXE']) + grxe_err = np.array(df0['GRXE_ERR']) + perc = np.percentile(grxe_err, grxe_err_cut, axis=0, keepdims=False) + print("{} {}: {}% cut of GRXE ERR: {:.2f} mCrab".format(key,enkey,grxe_err_cut,perc)) + df0 = df.query('DS9X == {} & DS9Y == {} & GRXE_ERR < {}'.format(ds9i,ds9j,perc)) + if (len(df0) <= nscw_min): + continue + + obsid = np.array(df0['OBSID']) grxe = np.array(df0['GRXE']) + grxe_err = np.array(df0['GRXE_ERR']) sg_mean, sg_med, sg_std = sigma_clipped_stats(grxe, sigma=sigma, maxiters=maxiters) filtered_data = sigma_clip(grxe, sigma=sigma, maxiters=maxiters, return_bounds=True) @@ -115,7 +131,8 @@ for i in range(sx): # final error on flux measurement ~RMS/sqrt(n) sg_sem = sem(grxe[filtered_grxe.mask==False]) - + obsid=obsid[filtered_grxe.mask==False] + #print("Sigma clipping: mean {:.2f} med {:.2f} std {:.2f} sem {:.2f}".format(sg_mean, sg_med, sg_std, sg_sem)) #plt.hist(grxe, bins=n_bins, range=minmax) @@ -126,8 +143,15 @@ for i in range(sx): sem_map[j][i] = sg_sem sign_map[j][i] = sg_mean/sg_sem cnt_map[j][i] = df0.shape[0] + + dkey="{:04d}{:04d}".format(j,i) + obsid_map[dkey] = obsid + grxe_map[dkey] = grxe + grxe_err_map[dkey] = grxe_err + +""" Filter by error map """ # Calculate the percentiles across the x and y dimension perc = np.percentile(sem_map, sem_cut, axis=(0, 1), keepdims=False) @@ -153,3 +177,52 @@ hdulist.writeto(mapsdir+fn.replace(".fits",".cnt.fits"),overwrite=True) hdulist[0].data=sign_map hdulist.writeto(mapsdir+fn.replace(".fits",".sign.fits"),overwrite=True) + +sys.exit() +print("Prepare data for fine map") + +obsid_list=[] +grxe_list=[] +grxe_err_list=[] +for i in range(sx): + for j in range(sy): + """ Use mask """ + if not (cnt_map[j][i]>0): + continue + world = w.wcs_pix2world([(i+1,j+1)], 1) + lon = world[0][0] + lat = world[0][1] + if(np.isnan(lon) or np.isnan(lat)): + continue + ds9i=i+1 + ds9j=j+1 + df0 = df.query('DS9X == {} & DS9Y == {}'.format(ds9i,ds9j)) + if (len(df0) <= nscw_min): + continue + + dkey="{:04d}{:04d}".format(j,i) + for scw in obsid_map[dkey]: + obsid_list.append(scw.decode("UTF-8")) + for grxe in grxe_map[dkey]: + grxe_list.append(grxe) + for grxe in grxe_err_map[dkey]: + grxe_err_list.append(grxe) + + +coldefs = fits.ColDefs([ + fits.Column(name='OBSID', format='11A', array=obsid_list), +]) + +fout = fn.replace(".fits",".grxe.fits") +hdu = fits.BinTableHDU.from_columns(coldefs, name='GRXE') +hdu.header['MISSION'] = ('INTEGRAL', '') +#hdu.header['TELESCOP'] = (outkey, '') +hdu.header['INSTITUT'] = ('IKI', 'Affiliation') +hdu.header['AUTHOR'] = ('Roman Krivonos', 'Responsible person') +hdu.header['EMAIL'] = ('krivonos@cosmos.ru', 'E-mail') +#hdu.add_checksum() +print(hdu.columns) +hdu.writeto(proddir+fout, overwrite=True) + +with fits.open(proddir+fout, mode='update') as hdus: + hdus[1].add_checksum() diff --git a/scripts/03_grxe_spec.py b/scripts/03_grxe_spec.py index 370d856..e8dc660 100755 --- a/scripts/03_grxe_spec.py +++ b/scripts/03_grxe_spec.py @@ -88,7 +88,7 @@ for skey in skeys: grxe_err = np.array(df['GRXE_ERR']) perc = np.percentile(grxe_err, grxe_err_cut, axis=0, keepdims=False) - print("{} {}: {}% cut of GRXE ERR: {:.2f} mCrab".format(skey,enkey,sem_cut,perc)) + print("{} {}: {}% cut of GRXE ERR: {:.2f} mCrab".format(skey,enkey,grxe_err_cut,perc)) idx=np.where(grxe_err < perc) grxe=grxe[idx] grxe_err=grxe_err[idx] @@ -105,9 +105,9 @@ for skey in skeys: print("{}: mean {:.2f} med {:.2f} std {:.2f} sem {:.2f} N={}".format(enkey, sg_mean, sg_med, sg_std, sg_sem, len(grxe[filtered_grxe.mask==False]))) #sg_sem*=1.5 - #if(sg_mean<0.0): - # sg_mean=1e-6 - # sg_sem*=2 + if(sg_mean<0.0): + sg_mean=1e-6 + #sg_sem*=2 ebands0[enkey]=[sg_mean,sg_sem]