SUBROUTINE projct(ear, ne, param, photar, photer) IMPLICIT NONE INCLUDE '../../inc/xspec.inc' INTEGER ne REAL ear(0:*), param(3), photar(ne), photer(ne) c Subroutine to do the 3-D to 2-D projection for a set of models in c shells onto a set of observations in annuli. The shells are assumed c to be prolate and to match the annuli. The major and minor axes for c the outer ellipse of each annulus are read from the XFLT### keywords c in the input files along with position angle and a set of start and c stop angles for the sections of the annulus included. The three c parameters are the semi-major and semi-minor axes for the inner ellipse c of the innermost annulus. This enables a hole to be left in the center c of the distribution (since shells only contribute to annuli outside c themselves). Obviously, these three parameters should not be left free ! c c 1/29/03 Modified to allow multiple observations each with their own set c of annuli. This requires that each observation have the same c number of annuli. If the observation 1 files are o1a1, o1a2, o1a3 c and the observation 2 files are o2a1, o2a2, o2a3 then they should c be read in using c XSPEC> data 1:1 o1a1 1:2 o2a1 2:3 o1a2 2:4 o2a2 3:5 o1a3 3:6 o2a3 c Arguments : c ear r i: Energy ranges c ne i i: 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 Pointers for the dynamic arrays INTEGER imajor, iminor, iorient, iprjmat INTEGER isects, iangbe, iangen INTEGER ismajor, isminor, isorient INTEGER issects, isangbe, isangen SAVE ismajor, isminor, isorient, issects, isangbe, isangen SAVE iprjmat c Local variables INTEGER nshells, nobs, i, j, ierr, snshells, snobs INTEGER maxsects SAVE snshells, snobs CHARACTER contxt*72 INTEGER dgndst, rgnenr, dgndtg, dgnflt EXTERNAL dgndst, rgnenr, dgndtg, dgnflt c The number of shells is the number of datagroups plus an inner "virtual" shell c to allow for a hole in the center if the first parameter (semi-major axis c of the hole) is greater than 0. The number of separate observations is the c the number of datasets divided by the number of datagroups. nshells = dgndtg() nobs = dgndst() / dgndtg() IF ( param(1) .GT. 0. ) nshells = nshells + 1 ierr = 0 c Check the number of filter keywords and calculate the max number of c sections maxsects = 1 DO i = 1, dgndst() IF ( dgnflt(i) .LT. 3 ) THEN WRITE(contxt, '(a,i2,a,i4)') & 'PROJCT: There are only ', dgnflt(i), & ' filter keywords for dataset ', i CALL xwrite(contxt, 10) WRITE(contxt, '(a)') & ' at least 3 are required.' CALL xwrite(contxt, 10) RETURN ENDIF maxsects = MAX(maxsects, (dgnflt(i)-3)/2) ENDDO c Grab the memory for the annulus descriptions CALL udmget(nshells*nobs, 7, imajor, ierr) contxt = 'Failed to grab memory for major axes' IF ( ierr .NE. 0 ) GOTO 999 CALL udmget(nshells*nobs, 7, iminor, ierr) contxt = 'Failed to grab memory for minor axes' IF ( ierr .NE. 0 ) GOTO 999 CALL udmget(nshells*nobs, 7, iorient, ierr) contxt = 'Failed to grab memory for orientations' IF ( ierr .NE. 0 ) GOTO 999 CALL udmget(nshells*nobs, 4, isects, ierr) contxt = 'Failed to grab memory for number of sections' IF ( ierr .NE. 0 ) GOTO 999 CALL udmget(nshells*maxsects*nobs, 7, iangbe, ierr) contxt = 'Failed to grab memory for start angle' IF ( ierr .NE. 0 ) GOTO 999 CALL udmget(nshells*maxsects*nobs, 7, iangen, ierr) contxt = 'Failed to grab memory for end angle' IF ( ierr .NE. 0 ) GOTO 999 c Grab the memory for the saved arrays and for the projection matrix CALL udmget(nshells*nobs, 7, ismajor, ierr) contxt = 'Failed to grab memory for saved major axes' IF ( ierr .NE. 0 ) GOTO 999 CALL udmget(nshells*nobs, 7, isminor, ierr) contxt = 'Failed to grab memory for saved minor axes' IF ( ierr .NE. 0 ) GOTO 999 CALL udmget(nshells*nobs, 7, isorient, ierr) contxt = 'Failed to grab memory for saved orientations' IF ( ierr .NE. 0 ) GOTO 999 CALL udmget(nshells*nobs, 4, issects, ierr) contxt = 'Failed to grab memory for saved number of sections' IF ( ierr .NE. 0 ) GOTO 999 CALL udmget(nshells*maxsects*nobs, 7, isangbe, ierr) contxt = 'Failed to grab memory for saved start angle' IF ( ierr .NE. 0 ) GOTO 999 CALL udmget(nshells*maxsects*nobs, 7, isangen, ierr) contxt = 'Failed to grab memory for saved end angle' IF ( ierr .NE. 0 ) GOTO 999 CALL udmget((nshells**2)*nobs, 7, iprjmat, ierr) contxt = 'Failed to grab memory for projection matrix' IF ( ierr .NE. 0 ) GOTO 999 c Run the main part of the routine. Split out to make the handling of c dynamic arrays easier CALL runprj(ne, ear, param, photar, photer, nshells, nobs, & MEMD(imajor), MEMD(iminor), & MEMD(iorient), maxsects, MEMI(isects), MEMD(iangbe), & MEMD(iangen), snshells, snobs, MEMD(ismajor), & MEMD(isminor), MEMD(isorient), MEMI(issects), & MEMD(isangbe), MEMD(isangen), MEMD(iprjmat)) c Release memory for arrays that don't have to be saved for the next c invocation. CALL udmfre(imajor, 7, ierr) contxt = 'Failed to free memory for major axes' IF ( ierr .NE. 0 ) GOTO 999 CALL udmfre(iminor, 7, ierr) contxt = 'Failed to free memory for minor axes' IF ( ierr .NE. 0 ) GOTO 999 CALL udmfre(iorient, 7, ierr) contxt = 'Failed to free memory for orientations' IF ( ierr .NE. 0 ) GOTO 999 CALL udmfre(isects, 4, ierr) contxt = 'Failed to free memory for number of sections' IF ( ierr .NE. 0 ) GOTO 999 CALL udmfre(iangbe, 7, ierr) contxt = 'Failed to free memory for start angle' IF ( ierr .NE. 0 ) GOTO 999 CALL udmfre(iangen, 7, ierr) contxt = 'Failed to free memory for end angle' IF ( ierr .NE. 0 ) GOTO 999 999 CONTINUE IF ( ierr .NE. 0 ) THEN CALL xwrite(contxt, 5) ENDIF RETURN END c ************************************************************************* SUBROUTINE runprj(ne, ear, param, photar, photer, nshells, nobs, & major, minor, orient, maxsects, sects, angbe, & angen, snshells, snobs, smajor, sminor, & sorient, ssects, sangbe, sangen, prjmat) IMPLICIT NONE INCLUDE '../../inc/xspec.inc' INTEGER ne, nshells, nobs, maxsects, snshells, snobs DOUBLE PRECISION major(nshells, nobs), minor(nshells, nobs) DOUBLE PRECISION orient(nshells, nobs) DOUBLE PRECISION angbe(maxsects, nshells, nobs) DOUBLE PRECISION angen(maxsects, nshells, nobs) DOUBLE PRECISION smajor(nshells, nobs), sminor(nshells, nobs) DOUBLE PRECISION sorient(nshells, nobs) DOUBLE PRECISION sangbe(maxsects, nshells, nobs) DOUBLE PRECISION sangen(maxsects, nshells, nobs) DOUBLE PRECISION prjmat(nshells, nshells, nobs) REAL ear(0:*), param(3), photar(ne), photer(ne) INTEGER sects(nshells, nobs) INTEGER ssects(nshells, nobs) INTEGER i, j, k, idt, ierr, nread INTEGER iphotmp CHARACTER fltstr*80, contxt*72 LOGICAL qcalc SAVE qcalc INTEGER dgnflt CHARACTER dgfilt*80 EXTERNAL dgnflt, dgfilt DATA qcalc / .TRUE. / nread = nshells IF ( param(1) .GT. 0. ) nread = nread - 1 c Load the annulus descriptions from the filter strings idt = 0 DO i = 1, nread DO j = 1, nobs idt = idt + 1 fltstr = dgfilt(1, idt) READ(fltstr, *) major(i,j) fltstr = dgfilt(2, idt) READ(fltstr, *) minor(i,j) fltstr = dgfilt(3, idt) READ(fltstr, *) orient(i,j) IF ( dgnflt(idt) .EQ. 3 ) THEN sects(i,j) = 1 angbe(1,i,j) = 0.d0 angen(1,i,j) = 360.d0 ELSE sects(i,j) = (dgnflt(idt)-3)/2 DO k = 1, sects(i,j) fltstr = dgfilt(3+2*k-1, idt) READ(fltstr, *) angbe(k,i,j) fltstr = dgfilt(3+2*k, idt) READ(fltstr, *) angen(k,i,j) ENDDO ENDIF ENDDO ENDDO c if required add the virtual shell from the three parameters IF ( param(1) .GT. 0. ) THEN DO j = 1, nobs major(nshells, j) = param(1) minor(nshells, j) = param(2) orient(nshells, j) = param(3) sects(nshells, j) = 1 angbe(1, nshells, j) = 0.d0 angen(1, nshells, j) = 360.d0 ENDDO ENDIF c Find out whether we need to calculate the projection matrix. We do c this if the number of shells have changed, if the annulus descriptions c have changed, or if this is first time through. IF ( .NOT. qcalc ) THEN IF ( nshells .NE. snshells ) qcalc = .TRUE. IF ( nobs .NE. snobs ) qcalc = .TRUE. DO i = 1, nshells DO j = 1, nobs IF ( major(i,j) .NE. smajor(i,j) .OR. & minor(i,j) .NE. sminor(i,j) .OR. & orient(i,j) .NE. sorient(i,j) .OR. & sects(i,j) .NE. ssects(i,j) & ) THEN qcalc = .TRUE. ELSE DO k = 1, sects(i,j) IF ( angbe(k,i,j) .NE. sangbe(k,i,j) .OR. & angen(k,i,j) .NE. sangen(k,i,j) ) THEN qcalc = .TRUE. ENDIF ENDDO ENDIF ENDDO ENDDO ENDIF IF ( qcalc ) THEN c Calculate the projection matrix DO j = 1, nobs CALL calprj(nshells, major(1,j), minor(1,j), orient(1,j), & maxsects, sects(1,j), angbe(1,1,j), & angen(1,1,j), prjmat(1,1,j)) ENDDO c Save the current values and reset the qcalc flag snshells = nshells snobs = nobs DO i = 1, nshells DO j = 1, nobs smajor(i,j) = major(i,j) sminor(i,j) = minor(i,j) sorient(i,j) = orient(i,j) ssects(i,j) = sects(i,j) DO k = 1, sects(i,j) sangbe(k,i,j) = angbe(k,i,j) sangen(k,i,j) = angen(k,i,j) ENDDO ENDDO ENDDO qcalc = .FALSE. ENDIF c Grab the memory for the temporary array required in the mixing of c the models. CALL udmget(ne, 6, iphotmp, ierr) contxt = 'Failed to grab memory for temporary model array' IF ( ierr .NE. 0 ) GOTO 999 c Calculate the projected models. CALL doproj(ne, ear, photar, photer, MEMR(iphotmp), nshells, & nobs, prjmat, major, minor) c Release the memory for the temporary array CALL udmfre(iphotmp, 6, ierr) contxt = 'Failed to free memory for temporary model array' IF ( ierr .NE. 0 ) GOTO 999 999 CONTINUE IF ( ierr .NE. 0 ) THEN CALL xwrite(contxt, 5) ENDIF RETURN END c ************************************************************************* SUBROUTINE calprj(nshells, major, minor, orient, maxsects, sects, & angbe, angen, prjmat) IMPLICIT NONE INTEGER nshells, maxsects DOUBLE PRECISION major(nshells), minor(nshells), orient(nshells) DOUBLE PRECISION angbe(maxsects, nshells) DOUBLE PRECISION angen(maxsects, nshells) DOUBLE PRECISION prjmat(nshells, nshells) INTEGER sects(nshells) c Subroutine to calculate the 3-D to 2-D projection matrix for prolate c ellipsoids onto elliptical shells. c Arguments : c nshells i: number of shells/annuli c major i: major axis length of annuli c minor i: minor axis length of annuli c orient i: orientation angle of annuli c maxsects i: max number of sections in any annulus c sects i: number of sections in each annulus c angbe i: start angle of annuli c angen i: end angle of annuli c prjmat r: projection matrix. note that element i,j c is the contribution of the ith ellipsoid to c the jth elliptical annulus. DOUBLE PRECISION majorm1 INTEGER idset, jdset, jdsetm1, i CHARACTER outstr*255 DOUBLE PRECISION prjint EXTERNAL prjint DO jdset = 1, nshells DO idset = 1, nshells prjmat(idset, jdset) = 0.d0 ENDDO ENDDO CALL xwrite('Projection matrix : ', 15) c Loop around target annuli DO jdset = 1, nshells c This is the target ellipse. Find the ellipse immediately interior to the c current one so that we can subtract off the contributions interior to the c defined annulus. jdsetm1 = 0 majorm1 = 0. DO i = 1, nshells IF ( major(i) .LT. major(jdset) ) THEN IF ( major(i) .GT. majorm1 ) THEN jdsetm1 = i majorm1 = major(i) ENDIF ENDIF ENDDO c and around contributing shells. DO idset = 1, nshells IF ( major(idset) .GE. major(jdset) ) THEN c Calculate the volume fraction of the current ellipsoid within the current c annulus. Note that we subtract off the volume within the inner edge c of the annulus - this is defined by the semi-major and semi-minor axes c and orientation of the next annulus but the angular sections of the c current annulus. prjmat(idset, jdset) & = prjint(major(idset), minor(idset), orient(idset), & major(jdset), minor(jdset), orient(jdset), & sects(jdset), angbe(1, jdset), & angen(1, jdset)) IF ( jdsetm1 .NE. 0 ) THEN prjmat(idset, jdset) = prjmat(idset, jdset) & - prjint(major(idset), minor(idset), orient(idset), & major(jdsetm1), minor(jdsetm1), & orient(jdsetm1), & sects(jdset), angbe(1, jdset), & angen(1, jdset)) ENDIF ENDIF ENDDO WRITE(outstr,'(25(1pg10.3))') & (prjmat(idset,jdset),idset=1,MIN(25,nshells)) CALL xwrite(outstr,15) ENDDO IF ( nshells .GT. 25 ) THEN CALL xwrite('Note: matrix rows truncated at 25 entries', 15) ENDIF RETURN END c ************************************************************************* FUNCTION prjint(shemaj, shemin, sheori, annmaj, annmin, annori, & annsec, annabe, annaen) INTEGER annsec DOUBLE PRECISION prjint, shemaj, shemin, sheori, annmaj DOUBLE PRECISION annmin, annori, annabe(annsec), annaen(annsec) c Calculate the fraction of the volume of the prolate ellipsoid with c semi-major axis shemaj, semi-minor axis shemin, and position angle c sheori which lies within the projected ellipse with semi-major axis c annmaj, semi-minor axis annmin, position angle annori, and within the c arc annabe to annaen as measured from the center. c Arguments : c prjint r: calculated fraction c shemaj i: shell semi-major axis c shemin i: shell semi-minor axis c sheori i: shell position angle c annmaj i: ellipse semi-major axis c annmin i: ellipse semi-minor axis c annori i: ellipse position angle c annsec i: number of ellipse sections c annabe i: ellipse start angle c annaen i: ellipse end angle c Integral to be calculated is c c V = 2/3 Int_(ph1-theta)^(ph2-theta) dx b^3 (1 + (b^2/a^2-1) cos^2 x)^-1 c c - 2/3 Int_(ph1-theta)^(ph2-theta) dx beta^3 c c (1 + (beta^2/alpha^2-1) cos^2 (x+theta-phi))^-3/2 c c ((b^2/beta^2)(1 + (beta^2/alpha^2-1) cos^2 (x+theta-phi)) - 1 c c - (b^2/a^2-1) cos^2 x)^3/2 (1 + (b^2/a^2-1) cos^2 x)^-1 c c where a,b,theta describe the ellipsoid and alpha,beta,phi,ph1,ph2 the c 2-D ellipse. INTEGER NSTEPS PARAMETER (NSTEPS=5000) DOUBLE PRECISION PI PARAMETER (PI=3.1415926535897932384626433d0) DOUBLE PRECISION b3, beta3, b2beta2, b2a2m1, beta2alpha2m1 DOUBLE PRECISION cos2x, cos2xtp, x, dx, f1, f2, angle1, angle2 DOUBLE PRECISION sumint, oridif INTEGER i, j c calculate constants that will be required in the integration. We c assume that the position angles input are in degrees so convert to c radians. prjint = 0.0d0 IF ( annmin .EQ. 0. .OR. shemaj .EQ. 0. OR. annmaj.EQ. 0. ) RETURN b3 = shemin**3 beta3 = annmin**3 b2beta2 = shemin**2/annmin**2 b2a2m1 = (shemin**2/shemaj**2 - 1) beta2alpha2m1 = (annmin**2/annmaj**2 - 1) oridif = (sheori - annori)*PI/180.d0 c Loop round sections DO j = 1, annsec angle1 = (annabe(j) - sheori)*PI/180.d0 angle2 = (annaen(j) - sheori)*PI/180.d0 c The integrand is smoothly varying so can use BFI and divide into NSTEP c steps from the start to end angles. dx = (angle2 - angle1) / NSTEPS sumint = 0.0d0 DO i = 1, NSTEPS x = dx * (i-0.5d0) + angle1 cos2x = COS(x)**2 cos2xtp = COS(x+oridif)**2 f1 = 1.0d0 + b2a2m1 * cos2x f2 = 1.0d0 + beta2alpha2m1 * cos2xtp IF ( (b2beta2*f2-f1) .NE. 0. ) THEN sumint = sumint + b3 / f1 - beta3 * f2**(-1.5d0) & * ( b2beta2 * f2 - f1 )**(1.5d0) / f1 ELSE sumint = sumint + b3 / f1 ENDIF ENDDO prjint = prjint + sumint * 2.0d0/3.0d0 * dx ENDDO RETURN END c ************************************************************************* SUBROUTINE doproj(ne, ear, photar, photer, photmp, nshells, & nobs, prjmat, major, minor) IMPLICIT NONE INTEGER ne, nshells, nobs REAL ear(0:ne), photar(ne), photer(ne), photmp(ne) DOUBLE PRECISION minor(nshells, nobs), major(nshells, nobs) DOUBLE PRECISION prjmat(nshells, nshells, nobs) c Mix up the models based on the projection matrix c Arguments : c ne i: total number of model bins c ear i: energy bins c photar i/r: input model values c photer i/r: input model errors c photmp w: workspace array c nshells i: number of shells/annuli (maybe the number c of datagroups + 1. The last shell is a virtual c shell at the center) c nobs i: number of separate observations c prjmat i: projection matrix. note that element i,j c is the contribution of the ith ellipsoid to c the jth elliptical annulus. c major i: semi-major axes of ellipses. c minor i: semi-minor axes of ellipses. INCLUDE '../../inc/xspec.inc' DOUBLE PRECISION PI PARAMETER(PI=3.1415926535897932384626433d0) DOUBLE PRECISION volfact, majorm1 INTEGER start, end, fstart, fend, photint INTEGER ie, i, idset, jdset, ishllm1, ierr INTEGER ishell, jshell, jobs INTEGER ine, iears, iphots INTEGER jne, jears, jphots CHARACTER contxt*72 INTEGER rgnenr, rgbear, rgdtse, dgndst EXTERNAL rgnenr, rgbear, rgdtse, dgndst ierr = 0 c Store the input model in the work array and initialize the output array DO ie = 1, ne photmp(ie) = photar(ie) photar(ie) = 0. ENDDO c Loop round the datasets DO jdset = 1, dgndst() c Set the number of energies and pointer in the ear and photar arrays c for this annulus jne = rgnenr(jdset) jears = rgbear(jdset) jphots = rgdtse(jdset)+1 c If this annulus has no energy bins then jump to the next annulus IF ( jne .EQ. 0 ) GOTO 200 c Set the observation number and shell number for this dataset jobs = mod(jdset-1, nobs) + 1 jshell = 1 + (jdset-1)/nobs c Now loop round the shells that are to be added together to give the c output spectrum DO ishell = 1, nshells idset = jobs + (ishell-1) * nobs c Find the ellipsoid interior to the current one so we can subtract off its c volume. ishllm1 = 0 majorm1 = 0.d0 DO i = 1, nshells IF ( major(i,jobs) .LT. major(ishell,jobs) ) THEN IF ( major(i,jobs) .GT. majorm1 ) THEN ishllm1 = i majorm1 = major(i,jobs) ENDIF ENDIF ENDDO c Volume factor is that of the current shell projected on the current c elliptical annulus divided by the total volume of the ellipsoidal c shell IF ( ishllm1 .NE. 0 ) THEN volfact = (prjmat(ishell, jshell, jobs) & -prjmat(ishllm1, jshell, jobs)) & /(4.0d0*PI/3.0d0) & /(major(ishell,jobs)*minor(ishell,jobs)**2 & -major(ishllm1,jobs)*minor(ishllm1,jobs)**2) ELSE volfact = prjmat(ishell, jshell, jobs) & /(4.0d0*PI/3.0d0) & /(major(ishell,jobs)*minor(ishell,jobs)**2) ENDIF IF ( volfact .GT. 0.0d0 .AND. rgnenr(idset) .NE. 0 ) THEN c If there is anything to multiply by then interpolate on the shell c spectrum to get the spectrum to add to the annulus. c Set the number of energies and pointer in the ear and photar arrays c for this shell ine = rgnenr(idset) iears = rgbear(idset) iphots = rgdtse(idset)+1 CALL udmget(jne, 4, start, ierr) contxt = 'Failed to get memory for start array' IF ( ierr .NE. 0 ) GOTO 999 CALL udmget(jne, 4, end, ierr) contxt = 'Failed to get memory for end array' IF ( ierr .NE. 0 ) GOTO 999 CALL udmget(jne, 6, fstart, ierr) contxt = 'Failed to get memory for fstart array' IF ( ierr .NE. 0 ) GOTO 999 CALL udmget(jne, 6, fend, ierr) contxt = 'Failed to get memory for fend array' IF ( ierr .NE. 0 ) GOTO 999 CALL udmget(jne, 6, photint, ierr) contxt = 'Failed to get memory for photint array' IF ( ierr .NE. 0 ) GOTO 999 CALL inibin(ine, ear(iears), jne, ear(jears), & MEMI(start), MEMI(end), MEMR(fstart), & MEMR(fend), 0.) CALL erebin(ine, photmp(iphots), jne, MEMI(start), & MEMI(end), MEMR(fstart), MEMR(fend), & MEMR(photint)) DO ie = 1, jne photar(jphots+ie-1) = photar(jphots+ie-1) & + MEMR(photint+ie-1) * SNGL(volfact) ENDDO CALL udmfre(start, 4, ierr) contxt = 'Failed to free memory from start array' IF ( ierr .NE. 0 ) GOTO 999 CALL udmfre(end, 4, ierr) contxt = 'Failed to free memory from end array' IF ( ierr .NE. 0 ) GOTO 999 CALL udmfre(fstart, 6, ierr) contxt = 'Failed to free memory from fstart array' IF ( ierr .NE. 0 ) GOTO 999 CALL udmfre(fend, 6, ierr) contxt = 'Failed to free memory from fend array' IF ( ierr .NE. 0 ) GOTO 999 CALL udmfre(photint, 6, ierr) contxt = 'Failed to free memory from photint array' IF ( ierr .NE. 0 ) GOTO 999 ENDIF ENDDO c Come from no energy bins in the current annulus 200 CONTINUE ENDDO 999 CONTINUE IF ( ierr .NE. 0 ) THEN CALL xwrite(contxt, 5) ENDIF RETURN END