SUBROUTINE xmmc(ear, netot, param, photar, photer) IMPLICIT NONE INCLUDE '../../inc/xspec.inc' INTEGER netot REAL ear(0:*), param(4), photar(netot), photer(netot) c Arguments : c ear r i: Energy ranges c netot i i: Total number of elements in photar array c param r i: Model parameters c photar r i/r: Model flux c photer r i/r: Model variance c Static arrays are c psfe r NPSFE INTEGER NPSFE PARAMETER(NPSFE=10) REAL psfe(NPSFE) c Dynamic arrays are c fact r nreg, nreg, npsfe, nobs c ireg i mxxsiz, mxysiz, nobs c surfbr r mxxsiz, mxysiz, nobs c image r imsiz(1), imsiz(2) c xsize i nobs c ysize i nobs c x0 r nobs c y0 r nobs c xmin i nobs c ymin i nobs c wmrbn i nobs c binsiz r nobs c instr i nobs c c where mxxsiz = max(xsize(nobs) and mxysiz = max(ysize(nobs)) INTEGER ifact, ine, iieoff, iireg, ixsize, iysize INTEGER ix0, iy0, ixmin, iymin, iwmrbn, ibinsiz, iwmap INTEGER iinstr, isurfbr, iimage DOUBLE PRECISION rapos, decpos, srapos, sdecpos INTEGER nreg, nobs, mxxsiz, mxysiz INTEGER snreg, snobs, snetot INTEGER status, i, imgsiz(2) CHARACTER contxt*72, imgfil*256, simgfil*256 LOGICAL qcentroid, qsurf, qforce SAVE srapos, sdecpos, snreg, snobs, snetot SAVE ifact, iireg, ixsize, iysize SAVE ix0, iy0, ixmin, iymin, iwmrbn, ibinsiz SAVE mxxsiz, mxysiz, psfe, iinstr, simgfil, isurfbr INTEGER dgndst, dgndtg EXTERNAL dgndst, dgndtg DATA snreg /0/ DATA psfe /0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 5.0, 6.0/ status = 0 c Check for any xset variables. If an image file is given then that will c be used for the surface brightness distribution. In this case also c return the image size. If no image file is given but a source center c RA and DEC are specified then these are used as the center of the surface c brightness distribution. If none of these are set the source center will c be taken as the wmap centroid. CALL xmmgps(imgfil, imgsiz, rapos, decpos, status) contxt = 'Failed to get image filename or position of source' IF ( status .NE. 0 ) GOTO 999 qcentroid = .FALSE. IF ( rapos .LT. 0. ) qcentroid = .TRUE. qsurf = .FALSE. IF ( imgfil .NE. ' ' ) qsurf = .TRUE. c Find the number of regions and observations - assume that each region is c in a different datagroup and that the number of datasets in each datagroup c gives the number of observations nreg = dgndtg() nobs = dgndst()/nreg c If things have changed then need to reallocate arrays and read information c from files qforce = .FALSE. IF ( netot .NE. snetot .OR. nreg .NE. snreg .OR. & nobs .NE. snobs .OR. ABS(rapos-srapos) .GT. 1.e-6 .OR. & ABS(decpos-sdecpos) .GT. 1.e-6 .OR. & imgfil .NE. simgfil ) THEN WRITE(contxt, '(i4, a, i4, a)') nreg, ' regions in ', nobs, & ' observations' CALL xwrite(contxt, 25) c Get the memory for the dynamic arrays whose size is the number of c observations CALL udmget(nobs, 4, ixsize, status) CALL udmget(nobs, 6, ix0, status) CALL udmget(nobs, 6, iy0, status) CALL udmget(nobs, 6, ixmin, status) CALL udmget(nobs, 6, iymin, status) CALL udmget(nobs, 4, iwmrbn, status) CALL udmget(nobs, 6, ibinsiz, status) CALL udmget(nobs, 4, iinstr, status) contxt = 'Failed to get memory for nobs arrays' IF ( status .NE. 0 ) GOTO 999 c Read information from the spectrum files to get the sizes of the wmaps, c and positions of the first wmap bin in detector coordinates. DO i = 0, nobs-1 CALL xmmgrg((i+1), nreg, nobs, MEMI(ixsize+i), & MEMI(iysize+i), MEMR(ixmin+i), MEMR(iymin+i), & MEMI(iwmrbn+i), MEMR(ibinsiz+i), & MEMI(iinstr+i), status) contxt = 'Failed to read spectrum files' IF ( status .NE. 0 ) GOTO 999 ENDDO c Find the max size of the WMAPS and number of energy bins mxxsiz = 0 mxysiz = 0 DO i = 0, nobs-1 mxxsiz = MAX(mxxsiz, MEMI(ixsize+i)) mxysiz = MAX(mxysiz, MEMI(iysize+i)) ENDDO c Allocate the memory for rest of the dynamic arrays CALL udmget((nreg**2)*NPSFE*nobs, 6, ifact, status) CALL udmget(mxxsiz*mxysiz, 4, iwmap, status) CALL udmget(mxxsiz*mxysiz*nobs, 4, iireg, status) CALL udmget(mxxsiz*mxysiz*nobs, 6, isurfbr, status) IF ( qsurf ) THEN CALL udmget(imgsiz(1)*imgsiz(2), 6, iimage, status) ENDIF contxt = 'Failed to get memory for second set of arrays' IF ( status .NE. 0 ) GOTO 999 c Load the ireg array and set the position of the source center in wmap bins DO i = 0, nobs-1 CALL xmmlrg((i+1), nreg, nobs, mxxsiz, mxysiz, & MEMI(ixsize+i), MEMI(iysize+i), & MEMR(ixmin+i), MEMR(iymin+i), MEMI(iwmrbn+i), & MEMR(ibinsiz+i), qcentroid, rapos, decpos, & MEMI(iwmap), MEMI(iireg+i*mxxsiz*mxysiz), & MEMR(ix0+i), MEMR(iy0+i), qsurf, & MEMR(isurfbr+i*mxxsiz*mxysiz), imgfil, & imgsiz(1), imgsiz(2), MEMR(iimage), status) ENDDO contxt = 'Failed to load ireg array' IF ( status .NE. 0 ) GOTO 999 CALL udmfre(iwmap, 4, status) IF ( qsurf ) CALL udmfre(iimage, 6, status) contxt = 'Failed to free memory for wmap and image work arrays' IF ( status .NE. 0 ) GOTO 999 snetot = netot snreg = nreg snobs = nobs srapos = rapos sdecpos = decpos simgfil = imgfil qforce = .TRUE. ENDIF c Run the mixing model CALL rnxmmc(netot, param, nobs, nreg, qforce, & MEMI(ixsize), MEMI(iysize), ear, & mxxsiz, mxysiz, MEMI(iireg), MEMR(ix0), MEMR(iy0), & MEMR(ixmin), MEMR(iymin), MEMI(iwmrbn), & MEMR(ibinsiz), MEMI(iinstr), NPSFE, psfe, qsurf, & MEMR(isurfbr), MEMR(ifact), photar, photer) 999 CONTINUE IF ( status .NE. 0 ) THEN CALL xwrite(contxt, 10) WRITE(contxt, '(a,i4)') 'XMMC : status = ', status CALL xwrite(contxt, 10) ENDIF RETURN END c ********************************************************************** SUBROUTINE rnxmmc(netot, param, nobs, nreg, qforce, xsize, & ysize, ear, mxxsiz, mxysiz, ireg, x0, & y0, xmin, ymin, wmrbn, binsiz, instr, npsfe, & psfe, qsurf, surfbr, fact, photar, photer) IMPLICIT NONE INCLUDE '../../inc/xspec.inc' INTEGER nobs, nreg, mxxsiz, mxysiz, netot, npsfe REAL ear(0:*), param(4), xmin(nobs), ymin(nobs), binsiz(nobs) REAL fact(nreg, nreg, npsfe, nobs), x0(nobs), y0(nobs) REAL photar(netot), photer(netot), psfe(npsfe) REAL surfbr(mxxsiz, mxysiz, nobs) INTEGER xsize(nobs), ysize(nobs) INTEGER ireg(mxxsiz, mxysiz, nobs), wmrbn(nobs) INTEGER instr(nobs) LOGICAL qforce, qsurf c Pointers for work arrays INTEGER ipsfarr, iwork, ivtmp c Local variables REAL alpha, beta, core INTEGER switch, status, iobs LOGICAL qnewsb, qerr CHARACTER outstr*256 SAVE alpha, beta, core, switch, qnewsb DATA qnewsb / .TRUE. / status = 0 c check whether the surface brightness parameter values have changed c and set them IF ( ABS(alpha-param(1)) .GT. 1.e-6 .OR. & ABS(beta -param(2)) .GT. 1.e-6 .OR. & ABS(core -param(3)) .GT. 1.e-6 .OR. & switch .NE. NINT(param(4)) .OR. & qforce ) qnewsb = .TRUE. alpha = param(1) beta = param(2) core = param(3) switch = NINT(param(4)) c If the surface brightness has changed then need to recalculate the c fact array for all the observations. IF ( qnewsb ) THEN CALL udmget(mxxsiz*mxysiz, 6, ipsfarr, status) DO iobs = 1, nobs WRITE(outstr, '(a,i4)') & 'Calculating mixing factors for observation ', iobs CALL xwrite(outstr, 25) CALL xmmfac( nreg, mxxsiz, mxysiz, xsize(iobs), & ysize(iobs), ireg(1,1,iobs), & x0(iobs), y0(iobs), xmin(iobs), ymin(iobs), & wmrbn(iobs), binsiz(iobs), instr(iobs), alpha, & beta, core, switch, qsurf, surfbr(1,1,iobs), & MEMR(ipsfarr), npsfe, psfe, fact(1,1,1,iobs)) ENDDO CALL udmfre(ipsfarr, 6, status) qnewsb = .FALSE. ENDIF c Mix the spectra CALL udmget(netot, 6, iwork, status) qerr = photer(netot/nobs/nreg/2) .GE. 0. IF ( qerr ) THEN CALL udmget(netot, 6, ivtmp, status) ENDIF CALL xmmmix(netot, nreg, nobs, ear, npsfe, psfe, & fact, MEMR(iwork), photar, qerr, MEMR(ivtmp), photer) CALL udmfre(iwork, 6, status) IF ( qerr ) CALL udmfre(ivtmp, 6, status) RETURN END c *********************************************************************** SUBROUTINE xmmfac(nreg, mxxsiz, mxysiz, xsize, ysize, & ireg, x0, y0, xmin, ymin, wmrbn, binsiz, instr, & alpha, beta, core, switch, qsurf, surfbr, & psfarr, npsfe, psfe, fact) IMPLICIT NONE INTEGER nreg, npsfe, mxxsiz, mxysiz, xsize, ysize REAL surfbr(mxxsiz, mxysiz), psfarr(mxxsiz, mxysiz) REAL fact(nreg, nreg, npsfe), psfe(npsfe) REAL xmin, ymin, x0, y0, alpha, beta, core, binsiz INTEGER ireg(mxxsiz, mxysiz) INTEGER wmrbn, instr, switch LOGICAL qsurf c nreg i i: number of regions c ireg i i: region containing this pixel (0=no region) c xmin, ymin r i: first wmap pixel in detector coordinates c wmrbn i i: wmap binning c binsiz r i: wmap bin size in arcmin c instr i i: instrument 1==MOS1, 2==MOS2, 3==PN c x0, y0 r i: center of source in wmap pixels c mxxsiz, mxysiz i i: size of the stored wmap array c xsize, ysize i i: actual size of the stored wmap c alpha r i: surface brightness parameter c beta r i: surface brightness parameter c core r i: surface brightness parameter (assumed arcmin) c switch i i: surface brightness parameter c qsurf l i: if true then surface brightness already calculated c surfbr i i: surface brightness c psfarr r w: workspace array for psf c npsfe i i: number of energies on which psf is calculated c psfe r i: energies on which psf is calculated c fact r r: weighting array c Routine to calculate weighting from source to target regions. Returns c fact(i,j,E) which is the weighting from region j to region i at energy c E. REAL corbin INTEGER ixs, iys, ixt, iyt, ie, iregs, iregt CHARACTER contxt*256 c Calculate the core radius in wmap bins corbin = core / binsiz c Calculate the normalized surface brightness array - normalization is c such that the surface brightness summed over each region is one. CALL xmmsb(mxxsiz, mxysiz, xsize, ysize, x0, y0, alpha, beta, & corbin, switch, nreg, ireg, qsurf, surfbr) c Initialize the fact array DO ie = 1, npsfe DO iregs = 1, nreg DO iregt = 1, nreg fact(iregt, iregs, ie) = 0.0 ENDDO ENDDO ENDDO c Loop over detector pixels DO iys = 1, ysize DO ixs = 1, xsize c If pixel is in a source region iregs = ireg(ixs, iys) IF ( iregs .EQ. 0 ) GOTO 100 c Loop over energies DO ie = 1, npsfe c Calculate PSF contribution to target detector pixels CALL xmmpsf(ixs, iys, xmin, ymin, wmrbn, binsiz, & mxxsiz, mxysiz, ireg, xsize, ysize, instr, & psfe(ie), psfarr) c Loop over target detector pixels DO iyt = 1, ysize DO ixt = 1, xsize c Accumulate in the fact array iregt = ireg(ixt, iyt) IF ( iregt .NE. 0 .AND. & psfarr(ixt, iyt) .GT. 0. ) THEN fact(iregt, iregs, ie) = fact(iregt, iregs, ie) & + surfbr(ixs,iys) * psfarr(ixt, iyt) ENDIF c End loop over target detector pixels ENDDO ENDDO c End loop over energies ENDDO 100 CONTINUE c End loop over source detector pixels ENDDO ENDDO DO ie = 1, npsfe WRITE(contxt,'(a,f7.3)') 'Energy = ', psfe(ie) CALL xwrite(contxt, 25) DO iregs = 1, nreg WRITE(contxt,*) (fact(iregt, iregs, ie), iregt=1,nreg) CALL xwrite(contxt, 25) ENDDO contxt = ' ' CALL xwrite(contxt, 25) ENDDO RETURN END c *********************************************************************** SUBROUTINE xmmgps(imgfil, imgsiz, rapos, decpos, status) IMPLICIT NONE DOUBLE PRECISION rapos, decpos INTEGER imgsiz(2), status CHARACTER imgfil*256 c Routine to get the source center RA and Dec (in degrees) as set by the c user using the xset command. If XMMPSF-RA has not been set then return c rapos=-999.0. c Arguments c imgfil c r: Image filename for surface brightness distribution c imgsiz(2) i r: Size of image c rapos d r: RA of source center in degrees c decpos d r: Dec of source center in degrees c status i r: 0==OK c 1==XMMPSF-RA read but no XMMPSF-DEC c 2==failed to read RA and Dec from string c -1==SLA error bad degrees c -2==SLA error bad arcmin c -3==SLA error bad arcsec c cfitsio error==failed to read image size DOUBLE PRECISION TPIINV PARAMETER (TPIINV=0.1591549430918953358D0) INTEGER iptr, iunit, findex CHARACTER instrg*256, cvalue*20, comment*60 CHARACTER fgmstr*256 EXTERNAL fgmstr status = 0 imgfil = ' ' rapos = -999.0 c Check for an image filename cvalue = 'XMMPSF-IMAGE' imgfil = fgmstr(cvalue) IF ( imgfil .NE. ' ' ) THEN CALL getlun(iunit) CALL xsftop(iunit, imgfil, 0, findex, status) CALL ftgkyj(iunit, 'NAXIS1', imgsiz(1), comment, status) CALL ftgkyj(iunit, 'NAXIS2', imgsiz(2), comment, status) CALL ftclos(iunit, status) CALL frelun(iunit) RETURN ENDIF c Check for an RA cvalue = 'XMMPSF-RA' instrg = ' ' instrg = fgmstr(cvalue) IF ( instrg .EQ. ' ' ) RETURN c If the string contains colons then the RA is in sexagesimal format c in hours, minutes, seconds. Otherwise assume it is in decimal degrees. IF ( index(instrg,':') .EQ. 0 ) THEN READ(instrg, *, iostat=status) rapos IF ( status .NE. 0 ) RETURN ELSE c Remove the : because SLA_DAFIN wants blank separated numbers DO iptr = 1, len(instrg) IF ( instrg(iptr:iptr) .EQ. ':' ) instrg(iptr:iptr) = ' ' ENDDO iptr = 1 CALL sla_dafin(instrg, iptr, rapos, status) IF ( status .NE. 0 ) RETURN c SLA_DAFIN converts to radians and assumes that the input numbers c are in degrees so adjust to make RA from h:m:s to degrees rapos = rapos * TPIINV * 15.d0 * 360.d0 ENDIF c Repeat for XMMPSF-DEC... cvalue = 'XMMPSF-DEC' instrg = ' ' instrg = fgmstr(cvalue) IF ( instrg .EQ. ' ' ) THEN status = 1 RETURN ENDIF IF ( index(instrg,':') .EQ. 0 ) THEN READ(instrg, *, iostat=status) decpos IF ( status .NE. 0 ) RETURN ELSE DO iptr = 1, len(instrg) IF ( instrg(iptr:iptr) .EQ. ':' ) instrg(iptr:iptr) = ' ' ENDDO iptr = 1 CALL sla_dafin(instrg, iptr, decpos, status) IF ( status .NE. 0 ) RETURN decpos = decpos * TPIINV * 360.d0 ENDIF RETURN END c ************************************************************************* SUBROUTINE xmmgrg(iobs, nreg, nobs, xsize, ysize, xmin, ymin, & wmrbn, binsiz, instr, status) IMPLICIT NONE INTEGER nreg, iobs, nobs, status REAL xmin, ymin, binsiz INTEGER xsize, ysize, wmrbn, instr c Routine to read the spectrum files to get the sizes of the wmaps, c and the positions of the first wmap bin in detector coordinates. c Arguments c iobs i i: The observation number to process c nreg i i: number of regions c nobs i i: number of observations c xsize i r: max wmap x size for this observation c ysize i r: max wmap y size for this observation c xmin r r: x first wmap bin in detector coordinates c ymin r r: y first wmap bin in detector coordinates c binsiz r r: The wmap binsize in arcminutes c instr i r: Instrument in use 1==MOS1, 2==MOS2, 3==PN c status i r: 0==OK REAL xhigh, yhigh, rval1, rval2 INTEGER ireg, idt, nx, ny, skybin INTEGER iunit, index CHARACTER infile*256, contxt*256, comment*72 CHARACTER cinstr*20 INTEGER lenact CHARACTER fgflnm*256 EXTERNAL lenact, fgflnm status = 0 c Loop round the regions reading the offsets and sizes of the wmaps xmin = 99999. ymin = 99999. xhigh = 0. yhigh = 0. DO ireg = 1, nreg idt = iobs + (ireg-1)*nobs CALL getlun(iunit) infile = fgflnm(idt) CALL xsftop(iunit, infile, 0, index, status) contxt = 'Failed to open '//infile(:lenact(infile)) IF ( status .NE. 0 ) GOTO 999 CALL ftgkyj(iunit, 'CDELT1P', wmrbn, comment, status) CALL ftgkyj(iunit, 'NAXIS1', nx, comment, status) CALL ftgkyj(iunit, 'NAXIS2', ny, comment, status) CALL ftgkye(iunit, 'CRVAL1P', rval1, comment, status) CALL ftgkye(iunit, 'CRVAL2P', rval2, comment, status) WRITE(contxt, '(a,i5)') & 'Failed to read WMAP keywords for dataset ', idt IF ( status .NE. 0 ) GOTO 999 xmin = MIN(xmin, rval1) ymin = MIN(ymin, rval2) xhigh = MAX(xhigh, rval1 + nx*wmrbn) yhigh = MAX(yhigh, rval2 + ny*wmrbn) CALL ftgkye(iunit, 'PIXSIZE', binsiz, comment, status) CALL ftgkyj(iunit, 'SKYBIN', skybin, comment, status) WRITE(contxt, '(a,i5)') & 'Failed to read pixel size keywords for dataset ', idt IF ( status .NE. 0 ) GOTO 999 CALL ftgkys(iunit, 'INSTRUME', cinstr, comment, status) WRITE(contxt, '(a,i5)') & 'Failed to read INSTRUME keyword for dataset ', idt IF ( status .NE. 0 ) GOTO 999 IF ( cinstr(1:5) .EQ. 'EMOS1' ) THEN instr = 1 ELSEIF ( cinstr(1:5) .EQ. 'EMOS2' ) THEN instr = 2 ELSEIF ( cinstr(1:3) .EQ. 'EPN' ) THEN instr = 3 ELSE contxt = 'Unknown instrument : '//cinstr CALL xwrite(contxt, 10) contxt = 'Defaulting to EMOS1' CALL xwrite(contxt, 10) instr = 1 ENDIF CALL ftclos(iunit, status) CALL frelun(iunit) ENDDO c Set the WMAP binsize in arcminutes binsiz = (60.0 * binsiz * wmrbn) / FLOAT(skybin) c Set the WMAP size xsize = (xhigh - xmin)/wmrbn + 1 ysize = (yhigh - ymin)/wmrbn + 1 WRITE(contxt,*) 'XMMGRG: ', xmin, ymin, xsize, ysize, binsiz CALL xwrite(contxt, 25) 999 CONTINUE IF ( status .NE. 0 ) THEN CALL xwrite(contxt, 10) WRITE(contxt, '(a, i4)') 'XMMGRG : status = ', status CALL xwrite(contxt, 10) ENDIF RETURN END c ************************************************************************* SUBROUTINE xmmlrg(iobs, nreg, nobs, mxxsiz, mxysiz, xsize, ysize, & xmin, ymin, wmrbn, binsiz, qcentroid, rapos, & decpos, wmap, ireg, x0, y0, qsurf, surfbr, & imgfil, imgsz1, imgsz2, image, status) IMPLICIT NONE INTEGER mxxsiz, mxysiz, imgsz1, imgsz2 REAL surfbr(mxxsiz, mxysiz), image(imgsz1, imgsz2) INTEGER ireg(mxxsiz, mxysiz), iobs, nreg, nobs, wmrbn, status INTEGER wmap(mxxsiz, mxysiz), xsize, ysize DOUBLE PRECISION rapos, decpos REAL xmin, ymin, binsiz, x0, y0 LOGICAL qcentroid, qsurf CHARACTER imgfil*(*) c Routine to load the IREG array based on the WMAP contents and set the c position of the source center in wmap coordinates. c Arguments : c iobs i i: Observation c nreg i i: Number of regions c nobs i i: Number of observations c mxxsiz i i: X size of wmap array c mxysiz i i: Y size of wmap array c xsize i i: Actual X size of wmap for this observation c ysize i i: Actual Y size of wmap for this observation c wmrbn i i: WMAP binning c binsiz r i: WMAP binsize in arcmin c xmin r i: X detector coordinate of first WMAP bin c ymin r i: Y detector coordinate of first WMAP bin c qcentroid l i: true if source center to come from centroiding c rapos d i: RA set by user c decpos d i: Dec set by user c wmap i w: Workspace array to hold WMAP values c ireg i r: region of given WMAP bin c x0 r r: x source center in wmap bins c y0 r r: y source center in wmap bins c qsurf l i: true if surface brightness to be read from image c surfbr r r: surface brightness map from image c imgfil c i: filename for surface brightness image c imgsz1 i i: image workspace x size c imgsz2 i i: image workspace y size c image r w: workspace array for image c status i r: 0==OK DOUBLE PRECISION PI PARAMETER (PI=3.1415926535897932384626433d0) DOUBLE PRECISION rfxval, rfyval, rfxpix, rfypix, rfxdlt, rfydlt DOUBLE PRECISION xpos, ypos, xpxli, ypxli, xpxlo, ypxlo, xh, yh DOUBLE PRECISION xrval, yrval, xrpix, yrpix, xinc, yinc, rot DOUBLE PRECISION wfxf0, wfyf0, wfxh0, wfyh0, xsign, ysign, wftheta DOUBLE PRECISION det2sky(6,3) REAL rval1, rval2, xcntrd, ycntrd, srfsum, wmapsum INTEGER ir, idt, iunit, index, i, j, instr INTEGER xoff, yoff, nx, ny, ipxl, jpxl CHARACTER infile*256, contxt*256, comment*72, coordtype*10 CHARACTER cinstr*20 LOGICAL qanyf INTEGER lenact CHARACTER fgflnm*256 EXTERNAL lenact, fgflnm SAVE det2sky DATA det2sky / & 2.416711943880075D+04, 2.529931630805351D+04, & 7.150000000000000D+01, -2.425000000000000D+02, & -1.0D0, 1.0D0, & 2.414449223345453D+04, 2.527086327870846D+04, & 6.950000000000000D+01, -1.260000000000000D+02, & -1.0D0, 1.0D0, & 2.547981387888674D+04, 2.313843360628563D+04, & -2.201500000000000D+03, -1.10100000000000D+03, & -1.0D0, 1.0D0/ status = 0 c Initialize ireg array DO j = 1, mxysiz DO i = 1, mxxsiz ireg(i,j) = 0 ENDDO ENDDO xcntrd = 0.0 ycntrd = 0.0 wmapsum = 0.0 c Loop round the regions DO ir = 1, nreg idt = iobs + (ir-1)*nobs CALL getlun(iunit) infile = fgflnm(idt) CALL xsftop(iunit, infile, 0, index, status) contxt = 'Failed to open '//infile(:lenact(infile)) IF ( status .NE. 0 ) GOTO 999 CALL ftgkyj(iunit, 'NAXIS1', nx, comment, status) CALL ftgkyj(iunit, 'NAXIS2', ny, comment, status) CALL ftgkye(iunit, 'CRVAL1P', rval1, comment, status) CALL ftgkye(iunit, 'CRVAL2P', rval2, comment, status) WRITE(contxt, '(a,i5)') & 'Failed to read WMAP keywords for dataset ', idt IF ( status .NE. 0 ) GOTO 999 c Get the reference WCS keywords which describe the sky coordinates CALL ftgkyd(iunit, 'REFXCRVL', rfxval, comment, status) CALL ftgkyd(iunit, 'REFYCRVL', rfyval, comment, status) CALL ftgkyd(iunit, 'REFXCRPX', rfxpix, comment, status) CALL ftgkyd(iunit, 'REFYCRPX', rfypix, comment, status) CALL ftgkyd(iunit, 'REFXCDLT', rfxdlt, comment, status) CALL ftgkyd(iunit, 'REFYCDLT', rfydlt, comment, status) WRITE(contxt, '(a,i5)') & 'Failed to read REF pointing keywords for dataset ', idt IF ( status .NE. 0 ) GOTO 999 c Try to read the EXT* keywords. If not load the values from the data c statement and read the PA_PNT to get the roll angle CALL ftgkyd(iunit, 'EXTXF0', wfxf0, comment, status) IF ( status .EQ. 0 ) THEN CALL ftgkyd(iunit, 'EXTYF0', wfyf0, comment, status) CALL ftgkyd(iunit, 'EXTXH0', wfxh0, comment, status) CALL ftgkyd(iunit, 'EXTYH0', wfyh0, comment, status) CALL ftgkyd(iunit, 'EXTXSIGN', xsign, comment, status) CALL ftgkyd(iunit, 'EXTYSIGN', ysign, comment, status) CALL ftgkyd(iunit, 'EXTTHETA', wftheta, comment, status) WRITE(contxt, '(a,i5)') & 'Failed to read EXT* keywords for dataset ', idt IF ( status .NE. 0 ) GOTO 999 ELSE status = 0 CALL ftgkyd(iunit, 'PA_PNT', wftheta, comment, status) CALL ftgkys(iunit, 'INSTRUME', cinstr, comment, status) WRITE(contxt, '(a,i5)') & 'Failed to read PA_PNT or INSTRUME keywords for dataset ', idt IF ( status .NE. 0 ) GOTO 999 c Check which instrument and set the angle - note each instrument has a c different orientation wrt the satellite roll. IF ( cinstr .EQ. 'EMOS1' ) THEN instr = 1 wftheta = -wftheta*PI/180.d0 ELSEIF ( cinstr .EQ. 'EMOS2' ) THEN instr = 2 wftheta = (-wftheta-90.d0)*PI/180.d0 ELSEIF ( cinstr .EQ. 'EPN' ) THEN instr = 3 wftheta = (-wftheta+90.d0)*PI/180.d0 ELSE comment = 'Invalid value of INSTRUME : '//cinstr status = -1 GOTO 999 ENDIF wfxf0 = det2sky(1,instr) wfyf0 = det2sky(2,instr) wfxh0 = det2sky(3,instr) wfyh0 = det2sky(4,instr) xsign = det2sky(5,instr) ysign = det2sky(6,instr) ENDIF CALL ftg2dj(iunit, 0, 0, mxxsiz, nx, ny, wmap, qanyf, status) WRITE(contxt, '(a,i5)') & 'Failed to read WMAP for dataset ', idt IF ( status .NE. 0 ) GOTO 999 CALL ftclos(iunit, status) CALL frelun(iunit) c Calculate the offset of this wmap into the ireg array xoff = (rval1 - xmin)/wmrbn yoff = (rval2 - ymin)/wmrbn WRITE(contxt,*) 'XMMLRG : ', ir, xoff, yoff CALL xwrite(contxt, 25) c Loop through wmap setting the region in the ireg array and the centroid c sums in case we need them DO j = 1, ny DO i = 1, nx IF ( wmap(i,j) .GE. 0 ) THEN IF ( ireg(i+xoff, j+yoff) .GT. 0 ) THEN WRITE(contxt, '(a,i5,i5,a,i4,i4)') & 'Region overlap at ', i, j, ' between ', & ireg(i+xoff, j+yoff), ir CALL xwrite(contxt, 25) ENDIF ireg(i+xoff, j+yoff) = ir xcntrd = xcntrd + wmap(i,j) * (i+xoff) ycntrd = ycntrd + wmap(i,j) * (j+yoff) wmapsum = wmapsum + wmap(i,j) ENDIF ENDDO ENDDO ENDDO IF ( wmapsum .NE. 0. ) THEN xcntrd = xcntrd / wmapsum ycntrd = ycntrd / wmapsum ENDIF c If we are getting the surface brightness map from an image then do that IF ( qsurf ) THEN c Read the image and the WCS data CALL getlun(iunit) CALL xsftop(iunit, imgfil, 0, index, status) contxt = 'Failed to open '//imgfil(:lenact(imgfil)) IF ( status .NE. 0 ) GOTO 999 CALL ftg2de(iunit, 0, 0, imgsz1, imgsz1, imgsz2, image, qanyf, & status) contxt = 'Failed to read image' IF ( status .NE. 0 ) GOTO 999 CALL ftgics(iunit, xrval, yrval, xrpix, yrpix, xinc, yinc, rot, & coordtype, status) contxt = 'Failed to read WCS data for image' IF ( status .NE. 0 ) GOTO 999 CALL ftclos(iunit, status) CALL frelun(iunit) c Load the surface brightness array. For each pixel in the array calculate c the sky pixel using the relations c xf = wfxf0 + xsign*(xh-wfxh0)*cos(wftheta) + ysign*(yh-wfyh0)*sin(wftheta) c yf = wfyf0 - xsign*(xh-wfxh0)*sin(wftheta) + ysign*(yh-wfyh0)*cos(wftheta) c Note that wftheta = -pa_pnt*pi/180 srfsum = 0.0 DO j = 1, ysize yh = (j-1)*wmrbn + ymin - wfyh0 DO i = 1, xsize xh = (i-1)*wmrbn + xmin - wfxh0 c Calculate the sky coordinate xpxli = wfxf0 + xsign*xh*cos(wftheta) & + ysign*yh*sin(wftheta) ypxli = wfyf0 - xsign*xh*sin(wftheta) & + ysign*yh*cos(wftheta) c Convert to RA and DEC using the REF* keyword information in the spectrum c header CALL ftwldp(xpxli, ypxli, rfxval, rfyval, rfxpix, rfypix, & rfxdlt, rfydlt, 0.0d0, '-TAN', xpos, ypos, & status) c Convert to input surface brightness image coordinates using the WCS c keywords from the image CALL ftxypx(xpos, ypos, xrval, yrval, xrpix, yrpix, & xinc, yinc, rot, coordtype, xpxlo, ypxlo, & status) ipxl = INT(xpxlo) jpxl = INT(ypxlo) cd cd WRITE(*,*) i, j, xpxli, ypxli, xpos, ypos, xpxlo, ypxlo cd IF ( ipxl .GT. 0 .AND. ipxl .LE. imgsz1 .AND. & jpxl .GT. 0 .AND. jpxl .LE. imgsz2 ) THEN surfbr(i,j) = image(ipxl, jpxl) srfsum = srfsum + surfbr(i,j) ELSE surfbr(i,j) = 0.0 ENDIF ENDDO ENDDO ELSE c Otherwise find the source center position in WMAP bin coordinates. If c qcentroid is set this is just the WMAP centroid. IF ( qcentroid ) THEN x0 = xcntrd y0 = ycntrd WRITE(contxt,*) 'XMMLRG : ', x0, y0 CALL xwrite(contxt, 25) ELSE c Otherwise convert the position given by the user into sky coordinates then c detector coordinates then wmap bin coordinates CALL ftxypx(rapos, decpos, rfxval, rfyval, rfxpix, rfypix, & rfxdlt, rfydlt, 0.0d0, '-TAN', xpos, ypos, & status) contxt = 'Failure in FTXYPX' IF ( status .NE. 0 ) GOTO 999 xh = wfxh0 + xsign*(xpos-wfxf0)*cos(wftheta) & - xsign*(ypos-wfyf0)*sin(wftheta) yh = wfyh0 + ysign*(xpos-wfxf0)*sin(wftheta) & + ysign*(ypos-wfyf0)*cos(wftheta) x0 = (SNGL(xh) - xmin)/wmrbn y0 = (SNGL(yh) - ymin)/wmrbn WRITE(contxt,*) 'XMMLRG : ', rapos, decpos, xpos, ypos, & wftheta, xh, yh, xmin, ymin, x0, y0 CALL xwrite(contxt, 25) ENDIF ENDIF 999 CONTINUE IF ( status .NE. 0 ) THEN CALL xwrite(contxt, 10) WRITE(contxt, '(a, i4)') 'XMMLRG : status = ', status CALL xwrite(contxt, 10) ENDIF RETURN END c ************************************************************************* SUBROUTINE xmmsb(mxxsiz, mxysiz, xsize, ysize, x0, y0, alpha, & beta, core, switch, nreg, ireg, qsurf, surfbr) IMPLICIT NONE INTEGER mxxsiz, mxysiz, xsize, ysize REAL surfbr(mxxsiz, mxysiz) REAL x0, y0, alpha, beta, core INTEGER ireg(mxxsiz, mxysiz), nreg, switch LOGICAL qsurf c Routine to calculate the surface brightness contribution to bin (ixs,iys) c for a source centered at (x0,y0). c Arguments : c mxxsiz i i: X size of wmap array c mxysiz i i: Y size of wmap array c xsize i i: X size of array to create c ysize i i: Y size of array to create c x0 r i: X position of source c y0 r i: Y position of source c alpha r i: surface brightness model parameter c beta r i: surface brightness model parameter c core r i: surface brightness model core radius (in wmap bins) c switch i i: surface brightness model switch c nreg i i: number of regions c ireg i i: map describing which pixels are in which region c qsurf l i: if true surface brightness has already been calculated c surfbr r i/r: surface brightness INTEGER ISUB REAL SUBINV PARAMETER (ISUB=1, SUBINV=1.0/ISUB) REAL x, y, dist2, snorm INTEGER ixs, iys, ix, iy, jx, jy, ir c If the surface brightness has not been calculated then do so IF ( .NOT. qsurf ) THEN c Loop round bins DO iys = 1, ysize y = iys - y0 DO ixs = 1, xsize x = ixs - x0 c Accumulate surface brightness integrating over the bin using bfi surfbr(ixs,iys) = 0.0 c Simple beta model IF ( switch .EQ. 0 ) THEn DO jx = 1, ISUB DO jy = 1, ISUB dist2 = (x-(jx-ISUB/2.-0.5)*SUBINV)**2 & +(y-(jy-ISUB/2.-0.5)*SUBINV)**2 surfbr(ixs,iys) = surfbr(ixs,iys) + & SUBINV*SUBINV* & (1+dist2/core**2)**(-3*beta+0.5) ENDDO ENDDO c Or the two slope model ELSEIF ( switch .EQ. 1 ) THEN DO jx = 1, ISUB DO jy = 1, ISUB dist2 = (x-(jx-ISUB/2.-0.5)*SUBINV)**2 & +(y-(jy-ISUB/2.-0.5)*SUBINV)**2 surfbr(ixs,iys) = surfbr(ixs,iys) + & SUBINV*SUBINV*(dist2**(-alpha)*exp(-dist2)+ & dist2**(-beta)*(1.0-exp(-dist2))) ENDDO ENDDO ENDIF ENDDO ENDDO ENDIF c Do the normalization. The sb for each region should normalize to one. DO ir = 1, nreg snorm = 0.0 DO iys = 1, ysize DO ixs = 1, xsize IF ( ireg(ixs,iys) .EQ. ir ) THEN snorm = snorm + surfbr(ixs,iys) ENDIF ENDDO ENDDO IF ( snorm .NE. 0. ) THEN DO iys = 1, ysize DO ixs = 1, xsize IF ( ireg(ixs,iys) .EQ. ir ) THEN surfbr(ixs,iys) = surfbr(ixs,iys) / snorm ENDIF ENDDO ENDDO ENDIF ENDDO RETURN END c ************************************************************************* SUBROUTINE xmmpsf(ixs, iys, xmin, ymin, wmrbn, binsiz, & mxxsiz, mxysiz, ireg, xsize, ysize, instr, & energy, psfarr) IMPLICIT NONE INTEGER mxxsiz, mxysiz, xsize, ysize REAL psfarr(mxxsiz, mxysiz), xmin, ymin, energy, binsiz INTEGER ireg(mxxsiz, mxysiz), ixs, iys, wmrbn, instr c Routine to return the PSF (in PSFARR) for a photon at position (IXS,IYS) c and energy ENERGY. c Arguments : c ixs i i: x source position (in wmap bins) c iys i i: y source position (in wmap bins) c xmin r i: x position of first wmap bin in detector coords c ymin r i: y position of first wmap bin in detector coords c wmrbn i i: wmap bin factor c binsiz r i: wmap binsize in arcmin c mxxsiz i i: x size of wmap array c mxysiz i i: y size of wmap array c ireg i i: array indicating which wmap bin lies in which region c xsize i i: x size of psf map to create c ysize i i: y size of psf map to create c instr i i: instrument 1==MOS1, 2==MOS2, 3=PN c energy r i: energy c psfarr r r: PSF c PSF data is from EPIC-MCT-TN-011 by Simona Ghizzardi (Oct 8, 2001) c The PSF is assumed to be a King profile c c PSF = A / (1 + (r/r_c)^2)^(-alpha) c c where A = (alpha-1)/r_c^2/PI to ensure the PSF integrates to 1 and c c r_c = r1(I) + r2(I)*E + r3(I)*Theta + r4(I)*E*Theta c alpha = a1(I) + a2(I)*E + a3(I)*Theta + a4(I)*E*Theta c c E is the energy (keV), Theta the off-axis angle (arcmin) and I is the c instrument. REAL PI PARAMETER (PI=3.1415962) REAL r1(3), r2(3), r3(3), r4(3) REAL a1(3), a2(3), a3(3), a4(3) REAL theta, rcore2, alpha, psfnrm, dist2 REAL psfsum INTEGER i, j, j2 SAVE r1, r2, r3, r4, a1, a2, a3, a4 DATA r1/5.074, 4.759, 6.636/ DATA r2/-0.236, -0.203, -0.305/ DATA r3/0.002, 0.014, -0.175/ DATA r4/-0.0180, -0.0229, -0.0067/ DATA a1/1.472, 1.411, 1.525/ DATA a2/-0.010, -0.005, -0.015/ DATA a3/-0.001, -0.001, -0.012/ DATA a4/-0.0016, -0.0002, -0.001/ c Calculate the off-axis angle of the source. Currently assume the optical c axis is at detector coordinates (0,0) theta = ((ixs * wmrbn + xmin)*binsiz/wmrbn)**2 + & ((iys * wmrbn + ymin)*binsiz/wmrbn)**2 theta = sqrt(theta) c Calculate core radius^2 and alpha for this energy and theta rcore2 = r1(instr) + r2(instr) * energy + r3(instr) * theta + & r4(instr) * energy * theta rcore2 = rcore2**2 alpha = a1(instr) + a2(instr) * energy + a3(instr) * theta + & a4(instr) * energy * theta c Convert the core radius from arcsec to wmap bins rcore2 = rcore2 / (binsiz*60.)**2 c Calculate the norm. Include the area of the wmap bin (which = 1). psfnrm = (alpha-1.0)/rcore2/PI c Loop round the output array psfsum = 0.0 DO j = 1, ysize j2 = (iys - j)**2 DO i = 1, xsize IF ( ireg(i, j) .NE. 0 ) THEN c Calculate the distance between the source and the target bin dist2 = (ixs - i)**2 + j2 c and the PSF contribution psfarr(i,j) = psfnrm * (1.0 + dist2/rcore2)**(-alpha) psfsum = psfsum + psfarr(i,j) ELSE psfarr(i,j) = 0.0 ENDIF ENDDO ENDDO cd cd WRITE(*,*) ixs, iys, energy, theta, SQRT(rcore2), alpha, psfsum cd RETURN END c ************************************************************************* SUBROUTINE xmmmix(netot, nreg, nobs, ear, npsfe, psfe, & fact, work, photar, qerr, vtmp, photer) IMPLICIT NONE INTEGER netot, nreg, nobs, npsfe REAL fact(nreg, nreg, npsfe, nobs), work(netot), photar(netot) REAL vtmp(netot), photer(netot), psfe(npsfe), ear(0:*) LOGICAL qerr c Routine to mix the model spectra based on the fact array c Arguments : c netot i i: Total number of energy bins c nreg i i: Number of regions c nobs i i: Number of observations c ear r i: Response energy ranges c npsfe i i: Number of energy bins in fact array c psfe r i: Energies of fact array bins c fact r i: Mixing array c work r w: Work array c photar r i/r: Model spectrum c qerr l i: True if model errors included c vtmp r w: Work array c photer r i/r: Model spectrum errors REAL energy, dele, weight INTEGER ie, idt, ir, iobs, jr, ie0, je0, ipsfe, iear0 INTEGER rgdtse, rgnenr, rgbear EXTERNAL rgdtse, rgnenr, rgbear c Store the input spectrum in the work array and initialize the output c array DO ie = 1, netot work(ie) = photar(ie) photar(ie) = 0.0 ENDDO c If necessary repeat for the variance IF ( qerr ) THEN DO ie = 1, netot vtmp(ie) = photer(ie) photer(ie) = 0.0 ENDDO ENDIF c Loop round observations DO iobs = 1, nobs c Loop round target regions DO ir = 1, nreg idt = iobs + (ir-1)*nobs c Set the offset into the photar array for the dataset ie0 = rgdtse(idt) c For this target dataset loop round the source regions DO jr = 1, nreg c Set the offset into the photar array for this source dataset je0 = rgdtse((jr-1)*nobs+iobs) c Loop round energies ipsfe = 0 DO ie = 1, rgnenr(idt) iear0 = rgbear(idt) + ie c Find interpolation information from fact energy bins to current energy energy = (ear(iear0-1)+ear(iear0))/2. IF ( ipsfe .LT. npsfe ) THEN IF ( energy .GT. psfe(ipsfe+1) ) THEN ipsfe = ipsfe + 1 IF ( ipsfe .LT. npsfe ) & dele = psfe(ipsfe+1) - psfe(ipsfe) ENDIF ENDIF IF ( ipsfe .EQ. 0 ) THEN weight = fact(ir, jr, 1, iobs) ELSEIF ( ipsfe .EQ. npsfe ) THEN weight = fact(ir, jr, npsfe, iobs) ELSE weight = ( (psfe(ipsfe+1)-energy) & *fact(ir, jr, ipsfe, iobs) + & (energy-psfe(ipsfe)) & *fact(ir, jr, ipsfe+1, iobs) ) / dele ENDIF photar(ie0+ie) = photar(ie0+ie) & + work(je0+ie) * weight IF ( qerr ) THEN photer(ie0+ie) = photer(ie0+ie) & + vtmp(je0+ie) * weight**2 ENDIF ENDDO ENDDO ENDDO ENDDO RETURN END