
      SUBROUTINE projct(ear, ne, param, photar, photer)

      INCLUDE '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 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
      INTEGER iphotmp

      SAVE ismajor, isminor, isorient, issects, isangbe, isangen
      SAVE iprjmat

c Local variables

      INTEGER nshells, i, j, ierr, neone, snshells
      INTEGER maxsects

      SAVE snshells

      CHARACTER contxt*72, fltstr*80

      LOGICAL qequal, qcalc

      SAVE qcalc

      INTEGER dgndst, rgnenr, dgnflt
      CHARACTER dgfilt*80
      EXTERNAL dgndst, rgnenr, dgfilt, dgnflt

      DATA qcalc / .TRUE. /

c The number of shells is the number of datasets 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

      nshells = dgndst()
      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, 6, imajor, ierr)
      contxt = 'Failed to grab memory for major axes'
      IF ( ierr .NE. 0 ) GOTO 999
      CALL udmget(nshells, 6, iminor, ierr)
      contxt = 'Failed to grab memory for minor axes'
      IF ( ierr .NE. 0 ) GOTO 999
      CALL udmget(nshells, 6, iorient, ierr)
      contxt = 'Failed to grab memory for orientations'
      IF ( ierr .NE. 0 ) GOTO 999
      CALL udmget(nshells, 4, isects, ierr)
      contxt = 'Failed to grab memory for number of sections'
      IF ( ierr .NE. 0 ) GOTO 999
      CALL udmget(nshells*maxsects, 6, iangbe, ierr)
      contxt = 'Failed to grab memory for start angle'
      IF ( ierr .NE. 0 ) GOTO 999
      CALL udmget(nshells*maxsects, 6, iangen, ierr)
      contxt = 'Failed to grab memory for end angle'
      IF ( ierr .NE. 0 ) GOTO 999

c Load the annulus descriptions from the filter strings

      DO i = 1, dgndst()

         fltstr = dgfilt(1, i)
         READ(fltstr, *) MEMR(imajor+i-1)

         fltstr = dgfilt(2, i)
         READ(fltstr, *) MEMR(iminor+i-1)

         fltstr = dgfilt(3, i)
         READ(fltstr, *) MEMR(iorient+i-1)

         IF ( dgnflt(i) .EQ. 3 ) THEN
            MEMI(isects+i-1) = 1
            MEMR(iangbe+(i-1)*maxsects) = 0.
            MEMR(iangen+(i-1)*maxsects) = 360.
         ELSE
            MEMI(isects+i-1) = (dgnflt(i)-3)/2
            DO j = 1, MEMI(isects+i-1)
               fltstr = dgfilt(3+2*j-1, i)
               READ(fltstr, *) MEMR(iangbe+(i-1)*maxsects+j-1)
               fltstr = dgfilt(3+2*j, i)
               READ(fltstr, *) MEMR(iangen+(i-1)*maxsects+j-1)
            ENDDO
         ENDIF

      ENDDO

c if required add the virtual shell from the three parameters

      IF ( param(1) .GT. 0. ) THEN
         MEMR(imajor+nshells-1) = param(1)
         MEMR(iminor+nshells-1) = param(2)
         MEMR(iorient+nshells-1) = param(3)
         MEMI(isects+nshells-1) = 1
         MEMR(iangbe+(nshells-1)*maxsects) = 0.
         MEMR(iangen+(nshells-1)*maxsects) = 360.
      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.

         DO i = 1, nshells
            IF ( MEMR(imajor+i-1) .NE. MEMR(ismajor+i-1) .OR.
     &           MEMR(iminor+i-1) .NE. MEMR(isminor+i-1) .OR.
     &           MEMR(iorient+i-1) .NE. MEMR(isorient+i-1) .OR.
     &           MEMI(isects+i-1) .NE. MEMI(issects+i-1)
     &         ) THEN
               qcalc = .TRUE.
            ELSE
               DO j = 1, MEMI(isects+i-1)
                  IF ( MEMR(iangbe+(i-1)*maxsects+j-1) .NE.
     &                 MEMR(isangbe+(i-1)*maxsects+j-1) .OR.
     &                 MEMR(iangen+(i-1)*maxsects+j-1) .NE.
     &                 MEMR(isangen+(i-1)*maxsects+j-1) ) THEN
                     qcalc = .TRUE.
                  ENDIF
               ENDDO
            ENDIF
         ENDDO

      ENDIF

      IF ( qcalc ) THEN

c Grab the memory for the saved arrays and for the projection matrix

         CALL udmget(nshells, 6, ismajor, ierr)
         contxt = 'Failed to grab memory for saved major axes'
         IF ( ierr .NE. 0 ) GOTO 999
         CALL udmget(nshells, 6, isminor, ierr)
         contxt = 'Failed to grab memory for saved minor axes'
         IF ( ierr .NE. 0 ) GOTO 999
         CALL udmget(nshells, 6, isorient, ierr)
         contxt = 'Failed to grab memory for saved orientations'
         IF ( ierr .NE. 0 ) GOTO 999
         CALL udmget(nshells, 4, issects, ierr)
         contxt = 'Failed to grab memory for saved number of sections'
         IF ( ierr .NE. 0 ) GOTO 999
         CALL udmget(nshells*maxsects, 6, isangbe, ierr)
         contxt = 'Failed to grab memory for saved start angle'
         IF ( ierr .NE. 0 ) GOTO 999
         CALL udmget(nshells*maxsects, 6, isangen, ierr)
         contxt = 'Failed to grab memory for saved end angle'
         IF ( ierr .NE. 0 ) GOTO 999

         CALL udmget(nshells**2, 6, iprjmat, ierr)
         contxt = 'Failed to grab memory for projection matrix'
         IF ( ierr .NE. 0 ) GOTO 999

c Calculate the projection matrix

         CALL calprj(nshells, MEMR(imajor), MEMR(iminor), MEMR(iorient),
     &               maxsects, MEMI(isects), MEMR(iangbe), MEMR(iangen),
     &               MEMR(iprjmat))

c Save the current values and reset the qcalc flag

         snshells = nshells
         DO i = 1, nshells
            MEMR(ismajor+i-1) = MEMR(imajor+i-1)
            MEMR(isminor+i-1) = MEMR(iminor+i-1)
            MEMR(isorient+i-1) = MEMR(iorient+i-1)
            MEMI(issects+i-1) = MEMI(isects+i-1)
            DO j = 1, MEMI(isects+i-1)
               MEMR(isangbe+(i-1)*maxsects+j-1) 
     &               = MEMR(iangbe+(i-1)*maxsects+j-1)
               MEMR(isangen+(i-1)*maxsects+j-1) 
     &               = MEMR(iangen+(i-1)*maxsects+j-1)
            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, 
     &            MEMR(iprjmat), MEMR(imajor), MEMR(iminor))

c Release the memory for the temporary arrays

      CALL udmfre(iphotmp, 6, ierr)
      contxt = 'Failed to free memory for temporary model array'
      IF ( ierr .NE. 0 ) GOTO 999
      CALL udmfre(imajor, 6, ierr)
      contxt = 'Failed to free memory for major axes'
      IF ( ierr .NE. 0 ) GOTO 999
      CALL udmfre(iminor, 6, ierr)
      contxt = 'Failed to free memory for minor axes'
      IF ( ierr .NE. 0 ) GOTO 999
      CALL udmfre(iorient, 6, 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, 6, ierr)
      contxt = 'Failed to free memory for start angle'
      IF ( ierr .NE. 0 ) GOTO 999
      CALL udmfre(iangen, 6, 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 calprj(nshells, major, minor, orient, maxsects, sects,
     &                  angbe, angen, prjmat)

      INTEGER nshells, maxsects
      REAL major(nshells), minor(nshells), orient(nshells)
      REAL angbe(maxsects, nshells), angen(maxsects, nshells)
      REAL 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.

      REAL majorm1
      INTEGER idset, jdset, jsect, jdsetm1, i

      CHARACTER outstr*255

      REAL prjint
      EXTERNAL prjint

      DO jdset = 1, nshells
         DO idset = 1, nshells
            prjmat(idset, jdset) = 0.
         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,nshells)
         CALL xwrite(outstr,15) 

      ENDDO

      RETURN
      END

c *************************************************************************

      FUNCTION prjint(shemaj, shemin, sheori, annmaj, annmin, annori,
     &                annsec, annabe, annaen)

      INTEGER annsec
      REAL prjint, shemaj, shemin, sheori, annmaj, annmin, annori
      REAL 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=1000)
      REAL PI
      PARAMETER (PI=3.1415962)

      REAL b3, beta3, b2beta2, b2a2m1, beta2alpha2m1, oridif
      REAL cos2x, cos2xtp, x, dx, f1, f2, angle1, angle2
      REAL sumint

      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.0

      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.

c Loop round sections

      DO j = 1, annsec

         angle1 = (annabe(j) - sheori)*PI/180.
         angle2 = (annaen(j) - sheori)*PI/180.

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.0
         DO i = 1, NSTEPS

            x = dx * (i-0.5) + angle1

            cos2x = COS(x)**2
            cos2xtp = COS(x+oridif)**2

            f1 = 1.0 + b2a2m1 * cos2x
            f2 = 1.0 + beta2alpha2m1 * cos2xtp

            sumint = sumint + b3 / f1 - beta3 * f2**(-1.5)
     &              * ( b2beta2 * f2 - f1 )**(1.5) / f1
         
         ENDDO

         prjint = prjint + sumint * 2./3. * dx

      ENDDO

      RETURN
      END

c *************************************************************************

      SUBROUTINE doproj(ne, ear, photar, photer, photmp, nshells, 
     &                  prjmat, major, minor)

      INTEGER ne, nshells

      REAL ear(0:ne), photar(ne), photer(ne), photmp(ne)
      REAL prjmat(nshells, nshells), major(nshells), minor(nshells)

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 datasets + 1. The last shell is a virtual
c                               shell at the center)
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 'xspec.inc'

      REAL PI
      PARAMETER(PI=3.1415962)

      REAL volfact, majorm1

      INTEGER start, end, fstart, fend, photint

      INTEGER ie, i, idset, jdset, idsetm1, ierr
      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 Now loop round the shells that are to be added together to give the 
c output spectrum

         DO idset = 1, dgndst()

c Find the ellipsoid interior to the current one so we can subtract off its
c volume.

            idsetm1 = 0
            majorm1 = 0.
            DO i = 1, nshells
               IF ( major(i) .LT. major(idset) ) THEN
                  IF ( major(i) .GT. majorm1 ) THEN
                     idsetm1 = i
                     majorm1 = major(i)
                  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 ( idsetm1 .NE. 0 ) THEN
               volfact = (prjmat(idset, jdset)-prjmat(idsetm1, jdset))
     &                   /(4.*PI/3.)/(major(idset)*minor(idset)**2
     &                               -major(idsetm1)*minor(idsetm1)**2)
            ELSE
               volfact = prjmat(idset, jdset)
     &                   /(4.*PI/3.)/(major(idset)*minor(idset)**2)
            ENDIF

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

            IF ( volfact .GT. 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.

               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) * 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

      ENDDO

 999  CONTINUE
      IF ( ierr .NE. 0 ) THEN
         CALL xwrite(contxt, 5)
      ENDIF

      RETURN
      END




